diff --git a/docs/UWGeodynamics/examples/1_23_08_CouplingModellingALEIB_TopographyRelaxation.ipynb b/docs/UWGeodynamics/examples/1_23_08_CouplingModellingALEIB_TopographyRelaxation.ipynb
new file mode 100644
index 00000000..920c28c3
--- /dev/null
+++ b/docs/UWGeodynamics/examples/1_23_08_CouplingModellingALEIB_TopographyRelaxation.ipynb
@@ -0,0 +1,197 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "9182138b-b954-4a5c-a6c6-35cb6c296923",
+   "metadata": {},
+   "source": [
+    "Topography Relaxation with surface processes\n",
+    "======\n",
+    "\n",
+    "This notebook models the topography relaxation with surface processes (Model TR-CM2) as described in Neng et al. (2026). It integrates Underworld2 with Badlands using the UWGeodynamic module within the ALE-IB scheme.\n",
+    "\n",
+    "**References**\n",
+    "\n",
+    "Lu, N., Moresi, L., Giordani, J., & Knight, B. (2026). A novel ALE scheme with the internal boundary for coupling tectonic and surface processes in geodynamic models. EGUsphere, 2026, 1-35."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "72d4536d-33cb-4c16-8812-645165b311f8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import underworld as uw\n",
+    "import math\n",
+    "from underworld import function as fn\n",
+    "import numpy as np\n",
+    "import os\n",
+    "\n",
+    "from underworld import UWGeodynamics as GEO\n",
+    "u = GEO.UnitRegistry\n",
+    "ndim = GEO.non_dimensionalise\n",
+    "dimen = GEO.dimensionalise\n",
+    "\n",
+    "comm = uw.mpi.comm\n",
+    "rank = uw.mpi.rank\n",
+    "size = uw.mpi.size\n",
+    "\n",
+    "GEO.rcParams[\"initial.nonlinear.tolerance\"] = 1e-3\n",
+    "GEO.rcParams['initial.nonlinear.max.iterations'] = 100\n",
+    "GEO.rcParams[\"nonlinear.tolerance\"] = 1e-3\n",
+    "GEO.rcParams['nonlinear.max.iterations'] = 100\n",
+    "GEO.rcParams[\"popcontrol.particles.per.cell.2D\"] = 20\n",
+    "GEO.rcParams[\"swarm.particles.per.cell.2D\"] = 20\n",
+    "GEO.rcParams[\"surface.pressure.normalization\"] = True\n",
+    "GEO.rcParams[\"pressure.smoothing\"] = True\n",
+    "GEO.rcParams[\"popcontrol.split.threshold\"] = 0.1\n",
+    "\n",
+    "half_rate = 1.0 * u.centimeter / u.year\n",
+    "model_length = 500. * u.kilometer\n",
+    "gravity = 9.81 * u.meter / u.second**2\n",
+    "bodyforce = 3300 * u.kilogram / u.metre**3 *gravity \n",
+    "\n",
+    "KL = model_length\n",
+    "Kt = KL / half_rate\n",
+    "KM = bodyforce * KL**2 * Kt**2\n",
+    "\n",
+    "GEO.scaling_coefficients[\"[length]\"] = KL\n",
+    "GEO.scaling_coefficients[\"[time]\"] = Kt\n",
+    "GEO.scaling_coefficients[\"[mass]\"]= KM"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d497849b-031f-4eee-af78-a5579245b274",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "longtest = False\n",
+    "dy = ndim(10.0 * u.kilometer)\n",
+    "max_time = 11.0*u.kiloyear\n",
+    "\n",
+    "if \"UW_LONGTEST\" in os.environ or longtest:\n",
+    "    dy = ndim(2.5 * u.kilometer)\n",
+    "    max_time = 301.0*u.kiloyear\n",
+    "\n",
+    "xmin, xmax = ndim(-160 * u.kilometer), ndim(160 * u.kilometer)\n",
+    "ymin, ymax = ndim(-160 * u.kilometer), ndim(40 * u.kilometer)\n",
+    "yint = 0.\n",
+    " \n",
+    "dy = ndim(2.5 * u.kilometer)\n",
+    "dx = dy\n",
+    "xRes,yRes = int(np.around((xmax-xmin)/dx)),int(np.around((ymax-ymin)/dy))\n",
+    "yResa,yResb =int(np.around((ymax-yint)/dy)),int(np.around((yint-ymin)/dy))\n",
+    "\n",
+    "Model = GEO.Model(elementRes=(xRes, yRes),\n",
+    "                  minCoord=(xmin,ymin),\n",
+    "                  maxCoord=(xmax, ymax),\n",
+    "                  gravity=(0.0, -gravity),\n",
+    "                  periodic=(True, False))\n",
+    "Model.outputDir= \"1_23_08_CouplingModellingALEIB_TopographyRelaxation_yres{:n}\".format(yRes)\n",
+    "Model.minStrainRate = 1e-18 / u.second\n",
+    "\n",
+    "wRatio = 1\n",
+    "D = np.abs(ymin)\n",
+    "Lambda = D/wRatio\n",
+    "k = 2.0 * np.pi / Lambda\n",
+    "mu0 = ndim(1e21  * u.pascal * u.second)\n",
+    "g = ndim(gravity)\n",
+    "rho0 = ndim(3300* u.kilogram / u.metre**3)\n",
+    "drho = rho0-0.\n",
+    "w_m = ndim(5*u.kilometer)\n",
+    "\n",
+    "tau0 = 2*k*mu0/drho/g\n",
+    "tau = (D*k+np.sinh(D*k)*np.cosh(D*k))/(np.sinh(D*k)**2)*tau0\n",
+    "\n",
+    "fn_coord = fn.input()\n",
+    "surf_fn = w_m * fn.math.cos(2.*np.pi*fn_coord[0]/Lambda) \n",
+    "\n",
+    "Model.inter_wall = Model._get_InternalwallSets(yint)\n",
+    "with Model.mesh.deform_mesh():\n",
+    "     Model.mesh.data[Model.inter_wall.data, 1] = surf_fn.evaluate(Model.inter_wall)[:,0]\n",
+    "Model._freeSurface_ALEIB = True \n",
+    "Model.freeSurface = True \n",
+    "Model._freeSurface.solve(0.)\n",
+    "\n",
+    "# recreate swarm as mesh deformed\n",
+    "import underworld as uw\n",
+    "from collections import OrderedDict \n",
+    "Model.swarm_variables = OrderedDict()\n",
+    "Model.swarm = uw.swarm.Swarm(mesh=Model.mesh, particleEscape=True)\n",
+    "Model.swarm.allow_parallel_nn = True \n",
+    "particlesPerCell = GEO.rcParams[\"swarm.particles.per.cell.2D\"]\n",
+    "Model._swarmLayout = uw.swarm.layouts.PerCellSpaceFillerLayout(swarm=Model.swarm,particlesPerCell=particlesPerCell)\n",
+    "Model.swarm.populate_using_layout(layout=Model._swarmLayout)\n",
+    "Model._initialize()\n",
+    "\n",
+    "materialAShape = fn_coord[1] > surf_fn\n",
+    "materialMShape = fn_coord[1] <= surf_fn\n",
+    "\n",
+    "materialA = Model.add_material(name=\"Air\", shape=materialAShape)\n",
+    "materialM = Model.add_material(name=\"Mantle\", shape=materialMShape) \n",
+    "\n",
+    "materialA.viscosity = 1e18 * u.pascal * u.second\n",
+    "materialM.viscosity = 1e21 * u.pascal * u.second\n",
+    "\n",
+    "materialA.density = 0.\n",
+    "materialM.density = 3300 * u.kilogram / u.metre**3\n",
+    "\n",
+    "sediment = Model.add_material(name=\"sediment\")\n",
+    "sediment.viscosity = 1e19 * u.pascal * u.second\n",
+    "sediment.density = 2700 * u.kilogram / u.metre**3 \n",
+    "\n",
+    "Model.set_velocityBCs(left=[0.,None],right=[0,None],bottom=[0.,0.], top=[None, 0.])\n",
+    "Model.init_model()\n",
+    "\n",
+    "Model.solver.set_inner_method(\"mumps\")\n",
+    "Model.solver.set_penalty(1e3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "72102da9-434f-4de3-b373-4c8db080694a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dt_set = 2.5*u.kiloyear\n",
+    "checkpoint_interval = 10.0*u.kiloyear\n",
+    "Model.surfaceProcesses = GEO.surfaceProcesses.Badlands(airIndex=[materialA.index],sedimentIndex=sediment.index,XML=\"resources/badlands_1en5.xml\", resolution=1.25 * u.kilometre, checkpoint_interval=dt_set,aspectRatio2d=0.25,surfElevation=surf_fn)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "daf7ba37-4acc-4054-9843-f9526ac32a9a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Model.run_for(max_time, checkpoint_interval=checkpoint_interval,dt=dt_set)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/docs/UWGeodynamics/examples/resources/badlands_1en5.xml b/docs/UWGeodynamics/examples/resources/badlands_1en5.xml
new file mode 100755
index 00000000..b9c3aad2
--- /dev/null
+++ b/docs/UWGeodynamics/examples/resources/badlands_1en5.xml
@@ -0,0 +1,111 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<badlands xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
+  <!-- Regular grid structure -->
+  <grid>
+    <!-- Optional parameter (integer) used to decrease TIN resolution.
+             The default value is set to 1. Increasing the factor
+             value will multiply the digital elevation model resolution
+             accordingly.  -->
+    <resfactor>1</resfactor>
+    <boundary>flat</boundary>
+    <!-- <boundary>fixed</boundary> -->
+    <!-- Underworld flag.  -->
+    <udw>1</udw>
+  </grid>
+  <!-- Sea-level structure -->
+  <sea>
+    <!-- Relative sea-level position [m] -->
+    <position>-0.</position>
+    <!-- Sea-level curve - (optional) -->
+    <!-- <curve>data/sea.csv</curve> -->
+    <!-- Limit flow network computation based on
+             water depth [m] -->
+    <!-- <limit>100.</limit> -->
+  </sea>
+  <!-- Simulation time structure -->
+  <time>
+    <!-- Simulation start time [a] -->
+    <start>0.</start>
+    <!-- Simulation end time [a] -->
+    <end>510000.</end>
+    <!-- Display interval [a] -->
+    <!-- This will be overridden by the linkage, so it doesn't matter -->
+    <display>10000.</display>
+  </time>
+  <strata>
+    <stratdx>1250.</stratdx>
+    <laytime>1250.</laytime>
+  </strata>
+  <!-- Precipitation structure -->
+  <precipitation>
+    <!-- Number of precipitation events -->
+    <climates>1</climates>
+    <!-- Precipitation definition -->
+    <rain>
+        <!-- Rain start time [a] -->
+        <rstart>0.</rstart>
+        <!-- Rain end time [a] -->
+        <rend>510000.</rend>
+        <!-- Rain computation time step [a] -->
+        <ortime>1250.</ortime>
+        <!-- Background precipitation value [m/a] -->
+        <rbgd>1.</rbgd>
+        <!-- Minimal precipitation value [m/a] -->
+        <rmin>0.1</rmin>
+        <!-- Maximal precipitation value [m/a] -->
+        <rmax>2.</rmax>
+        <!-- Wind velocity along X (W-E) direction [m/s] -->
+        <windx>2.</windx>
+        <!-- Wind velocity along Y (S-N) direction [m/s] -->
+        <windy>0.</windy>
+        <!-- Time conversion from cloud water to hydrometeors
+             range from 200 to 2000 [s]. Optional default is set
+             to 1000 s -->
+        <tauc>1000.</tauc>
+        <!-- Time for hydrometeor fallout range from 200 to 2000 [s].
+             Optional default is set to 1000 s -->
+        <tauf>1000.</tauf>
+        <!-- Moist stability frequency range from 0 to 0.01 [/s].
+             Optional default is set to 0.005 /s -->
+        <nm>0.005</nm>
+        <!-- Uplift sensitivity factor range from 0.001 to 0.02 [kg/m3].
+             Optional default is set to 0.005 kg/m3 -->
+        <cw>0.01</cw>
+        <!-- Depth of the moist layer range from 1000 to 5000 [m].
+             Optional default is set to 3000 m -->
+        <hw>5000.</hw>
+    </rain>
+  </precipitation>
+  <!-- Stream power law parameters:
+         The stream power law is a simplified form of the usual expression of
+         sediment transport by water flow, in which the transport rate is assumed
+         to be equal to the local carrying capacity, which is itself a function of
+         boundary shear stress. -->
+  <sp_law>
+    <!-- Values of m and n indicate how the incision rate scales
+             with bed shear stress for constant value of sediment flux
+             and sediment transport capacity.
+             Generally, m and n are both positive, and their ratio
+             (m/n) is considered to be close to 0.5 -->
+    <m>0.5</m>
+    <n>1.0</n>
+    <!-- The erodibility coefficient is scale-dependent and its value depend
+              on lithology and mean precipitation rate, channel width, flood
+              frequency, channel hydraulics. -->
+    <erodibility>1.e-5</erodibility>
+    <diffnb>10</diffnb>
+    <diffprop>0.05</diffprop>
+  </sp_law>
+  <!-- Linear slope diffusion parameters:
+         Parameterisation of the sediment transport includes the simple creep transport
+         law which states that transport rate depends linearly on topographic gradient. -->
+  <creep>
+    <!-- Surface diffusion coefficient [m2/a] -->
+    <caerial>0.1</caerial>
+    <!-- Marine diffusion coefficient [m2/a] -->
+    <cmarine>0.1</cmarine>
+    <criver>0.1</criver>
+  </creep>
+  <!-- Output folder path -->
+  <outfolder>outbdls</outfolder>
+</badlands>
diff --git a/src/underworld/UWGeodynamics/_freesurface.py b/src/underworld/UWGeodynamics/_freesurface.py
index 0bdc9584..be7b13e2 100644
--- a/src/underworld/UWGeodynamics/_freesurface.py
+++ b/src/underworld/UWGeodynamics/_freesurface.py
@@ -100,17 +100,34 @@ def _advect_surface(self, dt):
         uw.mpi.barrier()
         self.TField.syncronise()
 
-    def _update_mesh(self):
+    def _advect_surface_sp(self,dt):
+        surf_fn_badlands = self.model.surfaceProcesses.solve(dt)
+        if self.interface:
+            if self.model.mesh.dim == 2:
+                interpolate_x = self.model.mesh.data[self.interface.data,0]
+                interpolate_z = surf_fn_badlands(interpolate_x)
+                self.TField.data[self.interface.data, 0] = interpolate_z.copy()
+
+            if self.model.mesh.dim == 3:
+                interpolate_x = self.model.mesh.data[self.interface.data,0]
+                interpolate_x = self.model.mesh.data[self.interface.data,1]
+                interpolate_z = surf_fn_badlands((interpolate_x,interpolate_y))
+                self.TField.data[self.interface.data, 0] = interpolate_z.copy()
+        uw.mpi.barrier()
+        self.TField.syncronise()
 
+    def _update_mesh(self):
         with self.model.mesh.deform_mesh():
             # Last dimension is the vertical dimension
             self.model.mesh.data[:, -1] = self.TField.data[:, 0].copy()
 
     def solve(self, dtime):
         """ Advect free surface through dt and update the mesh """
-
         # First we advect the surface
-        self._advect_surface(dtime)
+        if self.model.surfaceProcesses:
+            self._advect_surface_sp(dtime)
+        else:
+            self._advect_surface(dtime)
         # Then we solve the system of linear equation
         self._solve_sle()
         # Finally we update the mesh
diff --git a/src/underworld/UWGeodynamics/surfaceProcesses.py b/src/underworld/UWGeodynamics/surfaceProcesses.py
index 16d95377..d2632d91 100644
--- a/src/underworld/UWGeodynamics/surfaceProcesses.py
+++ b/src/underworld/UWGeodynamics/surfaceProcesses.py
@@ -6,7 +6,7 @@
 import sys
 import math
 from scipy.ndimage.filters import gaussian_filter
-from scipy.interpolate import griddata, interp1d
+from scipy.interpolate import griddata, interp1d, CloughTocher2DInterpolator
 from underworld.scaling import non_dimensionalise as nd
 from underworld.scaling import dimensionalise
 from underworld.scaling import units as u
@@ -369,7 +369,7 @@ def solve(self, dt, sigma=0, uw_sample_style=0):
         self.time_years += dt_years
 
         # TODO: Improve the performance of this function
-        self._update_material_types()
+        surf_fn_badlands = self._update_material_types()
         comm.Barrier()
 
         if rank == 0 and self.verbose:
@@ -378,7 +378,7 @@ def solve(self, dt, sigma=0, uw_sample_style=0):
             print(purple + "Processing surface with Badlands...Done" + endcol)
             sys.stdout.flush()
 
-        return
+        return surf_fn_badlands if self.Model._freeSurface_ALEIB else None
 
     def _determine_particle_state_2D(self):
 
@@ -419,7 +419,7 @@ def _determine_particle_state_2D(self):
 
         flags = uw_surface[:, 1] < bdl_surface
 
-        return flags
+        return flags,f
 
     def _determine_particle_state(self):
         # Given Badlands' mesh, determine if each particle in 'volume' is above
@@ -463,15 +463,16 @@ def _determine_particle_state(self):
 
         # True for sediment, False for air
         flags = volume[:, 2] < interpolate_z
+        f = CloughTocher2DInterpolator((known_xy[:,0], known_xy[:,1]), known_z)
 
-        return flags
+        return flags,f
 
     def _update_material_types(self):
         # What do the materials (in air/sediment terms) look like now?
         if self.Model.mesh.dim == 3:
-            under_bd_surface = self._determine_particle_state()
+            under_bd_surface, surf_fn_badlands = self._determine_particle_state()
         if self.Model.mesh.dim == 2:
-            under_bd_surface = self._determine_particle_state_2D()
+            under_bd_surface,surf_fn_badlands = self._determine_particle_state_2D()
 
         # If any materials changed state, update the Underworld material types
         mi = self.Model.materialField.data
@@ -487,6 +488,7 @@ def _update_material_types(self):
             # if material is not air, and above surface, make it air
             eroded_mask = np.logical_and(~np.in1d(mi, air_material), ~under_bd_surface)
             mi[eroded_mask] = self.airIndex[0]
+        return surf_fn_badlands
 
     def _inject_badlands_displacement(self, time, dt, disp, sigma):
         """