diff --git a/NAMESPACE b/NAMESPACE
index 6f14efb..8e4893b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -22,6 +22,7 @@ export(GRIDobj2mat)
 export(GRIDobj2pm)
 export(GRIDobj2xyz)
 export(and)
+export(excesstopography)
 export(ezgetnal)
 export(fillsinks)
 export(flow_accumulation)
diff --git a/R/excesstopography.R b/R/excesstopography.R
new file mode 100644
index 0000000..3353ca6
--- /dev/null
+++ b/R/excesstopography.R
@@ -0,0 +1,69 @@
+#' Excess topography
+#'
+#' Reconstruct surface with threshold-slope surface
+#'
+#' The excess topography (Blöthe et al. 2015) is computed by solving
+#' an eikonal equation (Anand et al. 2023) constrained to lie below
+#' the original DEM. Where the slope of the DEM is greater than the
+#' threshold slope, the eikonal solver limits the output topography to
+#' that slope, but where the slope of the DEM is lower that the
+#' threshold slope, the output follows the DEM.
+#'
+#' The eikonal equation is solved using the fast sweeping method (Zhao
+#' 2004), which iterates over the DEM in alternating directions and
+#' updates the topography according to an upwind discretization of the
+#' gradient. To constrain the solution by the original DEM, the output
+#' topography is initiated with the DEM and only updates lower than
+#' the DEM are accepted.
+#'
+#' @param dem GRIDobj; Digital elevation model
+#' @param threshold_slope float | GRIDobj; The threshold
+#'     slope. Spatially variable slopes can be provided as a GRIDobj.
+#'
+#' @return GRIDobj; The reconstructed topography. The excess
+#'     topography is the difference between the original DEM and the
+#'     reconstructed one.
+#'
+#' @examples
+#'
+#' data(srtm_bigtujunga30m_utm11)
+#' DEM <- GRIDobj(srtm_bigtujunga30m_utm11)
+#' DEMext <- excesstopography(DEM, tan(20*pi/180))
+#'
+#' @export
+
+excesstopography <- function(dem, threshold_slope) {
+  # Input validation
+  dem <- processgrid(dem)
+  provided_go <- dem$provided_go
+  dem <- dem$r
+
+  # Preallocate output
+  ext <- dem
+
+  # Get DEM input
+  dem <- get_grid_data(dem)
+
+  # Pass grids to libtopotoolbox
+  if (is.numeric(threshold_slope)) {
+    threshold_grid <- as.single(rep(threshold_slope, length(dem$z)))
+  } else {
+    threshold_grid <- as.single(get_grid_data(processgrid(threshold_slope)$r)$z)
+  }
+  output <- single(length(dem$z))
+  result <- .C(
+    "wrap_excesstopography",
+    outputR = as.single(output),
+    demR = as.single(dem$z),
+    thresholdR = threshold_grid,
+    cellsizeR = as.single(dem$cellsize),
+    dimsR = as.integer(dem$dims),
+    NAOK = TRUE
+  )$outputR
+
+  terra::values(ext) <- result
+  if (provided_go) {
+    ext <- GRIDobj(ext)
+  }
+  return(ext)
+}
diff --git a/man/excesstopography.Rd b/man/excesstopography.Rd
new file mode 100644
index 0000000..6524d56
--- /dev/null
+++ b/man/excesstopography.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/excesstopography.R
+\name{excesstopography}
+\alias{excesstopography}
+\title{Excess topography}
+\usage{
+excesstopography(dem, threshold_slope)
+}
+\arguments{
+\item{dem}{GRIDobj; Digital elevation model}
+
+\item{threshold_slope}{float | GRIDobj; The threshold
+slope. Spatially variable slopes can be provided as a GRIDobj.}
+}
+\value{
+GRIDobj; The reconstructed topography. The excess
+    topography is the difference between the original DEM and the
+    reconstructed one.
+}
+\description{
+Reconstruct surface with threshold-slope surface
+}
+\details{
+The excess topography (Blöthe et al. 2015) is computed by solving
+an eikonal equation (Anand et al. 2023) constrained to lie below
+the original DEM. Where the slope of the DEM is greater than the
+threshold slope, the eikonal solver limits the output topography to
+that slope, but where the slope of the DEM is lower that the
+threshold slope, the output follows the DEM.
+
+The eikonal equation is solved using the fast sweeping method (Zhao
+2004), which iterates over the DEM in alternating directions and
+updates the topography according to an upwind discretization of the
+gradient. To constrain the solution by the original DEM, the output
+topography is initiated with the DEM and only updates lower than
+the DEM are accepted.
+}
+\examples{
+
+data(srtm_bigtujunga30m_utm11)
+DEM <- GRIDobj(srtm_bigtujunga30m_utm11)
+DEMext <- excesstopography(DEM, tan(20*pi/180))
+
+}
diff --git a/src/wrap_excesstopography.c b/src/wrap_excesstopography.c
new file mode 100644
index 0000000..d2ca63f
--- /dev/null
+++ b/src/wrap_excesstopography.c
@@ -0,0 +1,15 @@
+#include <R.h>
+#include <stddef.h>
+#include <stdint.h>
+
+#include "topotoolbox.h"
+#include "topotoolboxr.h"
+
+void wrap_excesstopography(float *excess, float *dem, float *threshold_slopes,
+                           float *cellsize, int *dimsR) {
+  // Transformation of integers and array allocation
+  ptrdiff_t dims[2] = {dimsR[0], dimsR[1]};
+
+  // Call libtopotoolbox
+  excesstopography_fsm2d(excess, dem, threshold_slopes, *cellsize, dims);
+}
diff --git a/tests/testthat/test-excesstopography.R b/tests/testthat/test-excesstopography.R
new file mode 100644
index 0000000..984c943
--- /dev/null
+++ b/tests/testthat/test-excesstopography.R
@@ -0,0 +1,12 @@
+test_that("Tests excesstopography", {
+  dem <- terra::rast(system.file("ex/elev.tif", package = "terra"))
+  dem <- terra::project(dem, "EPSG:32632", res = 90.0)
+  expect_silent(dem_ext <- excesstopography(GRIDobj(dem), tan(20 * pi / 180)))
+  expect_true(inherits(dem_ext, "GRIDobj"))
+  expect_true(all(terra::values(dem_ext$raster) <= terra::values(dem), na.rm=TRUE))
+
+  ts <- dem
+  terra::values(ts) <- runif(terra::size(ts))
+  expect_silent(dem_ext <- excesstopography(GRIDobj(dem), GRIDobj(ts)))
+  expect_true(all(terra::values(dem_ext$raster) <= terra::values(dem), na.rm=TRUE))
+})