From 09fad83d37e5b1225a74000b0d46f93eb89dd1c3 Mon Sep 17 00:00:00 2001
From: Sylwester Arabas <sylwester.arabas@agh.edu.pl>
Date: Tue, 30 Sep 2025 04:37:05 +0200
Subject: [PATCH 01/18] early draft of deterministic deposition removal

---
 PySDM/dynamics/__init__.py                    |  1 +
 PySDM/dynamics/deposition_removal.py          | 26 +++++++++
 PySDM/particulator.py                         | 34 +++++++++++
 docs/bibliography.json                        |  7 +++
 .../dynamics/test_deposition_removal.py       | 56 +++++++++++++++++++
 5 files changed, 124 insertions(+)
 create mode 100644 PySDM/dynamics/deposition_removal.py
 create mode 100644 tests/unit_tests/dynamics/test_deposition_removal.py

diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py
index d296bf4720..7a1fb1c354 100644
--- a/PySDM/dynamics/__init__.py
+++ b/PySDM/dynamics/__init__.py
@@ -17,3 +17,4 @@
 from PySDM.dynamics.relaxed_velocity import RelaxedVelocity
 from PySDM.dynamics.seeding import Seeding
 from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce
+from PySDM.dynamics.deposition_removal import DepositionRemoval
diff --git a/PySDM/dynamics/deposition_removal.py b/PySDM/dynamics/deposition_removal.py
new file mode 100644
index 0000000000..3b1a582734
--- /dev/null
+++ b/PySDM/dynamics/deposition_removal.py
@@ -0,0 +1,26 @@
+"""deposition removal logic for zero-dimensional environments"""
+
+from PySDM.dynamics.impl import register_dynamic
+import numpy as np
+
+
+@register_dynamic()
+class DepositionRemoval:
+    def __init__(self, *, all_or_nothing: bool):
+        """stochastic ("all or nothing") or deterministic (multiplicity altering) removal"""
+        self.all_or_nothing = all_or_nothing
+        self.length_scale = None
+
+    def register(self, builder):
+        builder.request_attribute("relative fall velocity")
+        assert builder.particulator.environment.mesh.n_dims == 0
+        self.particulator = builder.particulator
+        self.length_scale = np.cbrt(self.particulator.environment.mesh.dv)
+
+    def __call__(self):
+        """see, e.g., the naive scheme in Algorithm 1 in
+        [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)"""
+        self.particulator.deposition_removal(
+            all_or_nothing=self.all_or_nothing,
+            length_scale=self.length_scale,
+        )
diff --git a/PySDM/particulator.py b/PySDM/particulator.py
index 5919feb805..710f889d66 100644
--- a/PySDM/particulator.py
+++ b/PySDM/particulator.py
@@ -565,3 +565,37 @@ def thaw_instantaneous(self):
             cell=self.attributes["cell id"],
             temperature=self.environment["T"],
         )
+
+    def deposition_removal(self, *, all_or_nothing, length_scale):
+        # TODO: to be moved into backend (and attributes?)
+
+        prob_zero_events = self.formulae.trivia.poissonian_avoidance_function
+
+        def backend_deposition_removal_all_or_nothing(
+            *, relative_fall_velocity, length_scale, timestep, u01
+        ):
+            pass
+
+        def backend_deposition_removal_deterministic(
+            *, relative_fall_velocity, multiplicity, length_scale, timestep
+        ):
+            for i, velocity in enumerate(relative_fall_velocity):
+                removal_rate = velocity / length_scale
+                survive_prob = prob_zero_events(r=removal_rate, dt=timestep)
+                assert 0 <= survive_prob <= 1
+                multiplicity[i] *= survive_prob
+
+        if all_or_nothing:
+            backend_deposition_removal_all_or_nothing(
+                relative_fall_velocity=self.attributes["relative fall velocity"].data,
+                length_scale=length_scale,
+                timestep=self.dt,
+                u01=None,
+            )
+        else:
+            backend_deposition_removal_deterministic(
+                relative_fall_velocity=self.attributes["relative fall velocity"].data,
+                multiplicity=self.attributes["multiplicity"].data,
+                length_scale=length_scale,
+                timestep=self.dt,
+            )
diff --git a/docs/bibliography.json b/docs/bibliography.json
index deb122ac9e..0aad561a4c 100644
--- a/docs/bibliography.json
+++ b/docs/bibliography.json
@@ -973,5 +973,12 @@
     ],
     "label": "Barthazy and Schefold 2006 (Atmos. Res. 82)",
     "title": "Fall velocity of snowflakes of different riming degree and crystal types"
+  },
+  "https://doi.org/10.1016/j.jcp.2016.06.029": {
+    "usages": [
+      "PySDM/dynamics/deposition_removal.py"
+    ],
+    "label": "Curtis et al. 2016 (J. Comp. Phys. 322)",
+    "title": "Accelerated simulation of stochastic particle removal\nprocesses in particle-resolved aerosol models"
   }
 }
diff --git a/tests/unit_tests/dynamics/test_deposition_removal.py b/tests/unit_tests/dynamics/test_deposition_removal.py
new file mode 100644
index 0000000000..58d127c027
--- /dev/null
+++ b/tests/unit_tests/dynamics/test_deposition_removal.py
@@ -0,0 +1,56 @@
+from PySDM.dynamics import DepositionRemoval
+from PySDM.physics import si
+from PySDM.environments import Box
+from PySDM.builder import Builder
+from PySDM.backends import CPU
+from PySDM.products import ParticleConcentration, SuperDropletCountPerGridbox, Time
+from matplotlib import pyplot
+import pytest
+import numpy as np
+
+
+class TestDepositionRemoval:
+    @staticmethod
+    @pytest.mark.parametrize("all_or_nothing", (True, False))
+    def test_convergence_wrt_dt(all_or_nothing, plot=True):
+        # arrange
+        dt = 1 * si.s
+        dv = 666 * si.m**3
+        n_steps = 100
+        multiplicity = [1, 10, 100, 1000]
+        water_mass = [1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug]
+        backend_instance = CPU()
+
+        builder = Builder(
+            n_sd=len(multiplicity),
+            environment=Box(dv=dv, dt=dt),
+            backend=backend_instance,
+        )
+        builder.add_dynamic(DepositionRemoval(all_or_nothing=all_or_nothing))
+        particulator = builder.build(
+            attributes={
+                "multiplicity": np.asarray(multiplicity),
+                "signed water mass": np.asarray(water_mass),
+            },
+            products=(ParticleConcentration(), SuperDropletCountPerGridbox(), Time()),
+        )
+
+        # act
+        output = {name: [] for name in particulator.products}
+        for step in range(n_steps):
+            if step != 0:
+                particulator.run(steps=1)
+            for name, product in particulator.products.items():
+                output[name].append(product.get() + 0)
+
+        # plot
+        pyplot.title(f"{all_or_nothing=}")
+        pyplot.xlabel("time [s]")
+        pyplot.ylabel("particle concentration [m$^{-3}$]")
+        pyplot.plot(output["time"], output["particle concentration"])
+        if plot:
+            pyplot.show()
+        else:
+            pyplot.clf()
+
+        # assert

From 6f9a3045750a2672f15fc11c03c58a42574e68a8 Mon Sep 17 00:00:00 2001
From: Sylwester Arabas <sylwester.arabas@agh.edu.pl>
Date: Tue, 30 Sep 2025 07:08:02 +0200
Subject: [PATCH 02/18] renaming deposition -> sedimentation; plot=False not to
 hang CI

---
 PySDM/dynamics/__init__.py                           |  2 +-
 ...eposition_removal.py => sedimentation_removal.py} |  4 ++--
 PySDM/particulator.py                                | 10 +++++-----
 ...tion_removal.py => test_sedimentation_removal.py} | 12 +++++++-----
 4 files changed, 15 insertions(+), 13 deletions(-)
 rename PySDM/dynamics/{deposition_removal.py => sedimentation_removal.py} (92%)
 rename tests/unit_tests/dynamics/{test_deposition_removal.py => test_sedimentation_removal.py} (83%)

diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py
index 7a1fb1c354..45725a2dac 100644
--- a/PySDM/dynamics/__init__.py
+++ b/PySDM/dynamics/__init__.py
@@ -17,4 +17,4 @@
 from PySDM.dynamics.relaxed_velocity import RelaxedVelocity
 from PySDM.dynamics.seeding import Seeding
 from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce
-from PySDM.dynamics.deposition_removal import DepositionRemoval
+from PySDM.dynamics.sedimentation_removal import SedimentationRemoval
diff --git a/PySDM/dynamics/deposition_removal.py b/PySDM/dynamics/sedimentation_removal.py
similarity index 92%
rename from PySDM/dynamics/deposition_removal.py
rename to PySDM/dynamics/sedimentation_removal.py
index 3b1a582734..bdb3414f7e 100644
--- a/PySDM/dynamics/deposition_removal.py
+++ b/PySDM/dynamics/sedimentation_removal.py
@@ -5,7 +5,7 @@
 
 
 @register_dynamic()
-class DepositionRemoval:
+class SedimentationRemoval:
     def __init__(self, *, all_or_nothing: bool):
         """stochastic ("all or nothing") or deterministic (multiplicity altering) removal"""
         self.all_or_nothing = all_or_nothing
@@ -20,7 +20,7 @@ def register(self, builder):
     def __call__(self):
         """see, e.g., the naive scheme in Algorithm 1 in
         [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)"""
-        self.particulator.deposition_removal(
+        self.particulator.sedimentation_removal(
             all_or_nothing=self.all_or_nothing,
             length_scale=self.length_scale,
         )
diff --git a/PySDM/particulator.py b/PySDM/particulator.py
index 710f889d66..fbf89de93a 100644
--- a/PySDM/particulator.py
+++ b/PySDM/particulator.py
@@ -566,17 +566,17 @@ def thaw_instantaneous(self):
             temperature=self.environment["T"],
         )
 
-    def deposition_removal(self, *, all_or_nothing, length_scale):
+    def sedimentation_removal(self, *, all_or_nothing, length_scale):
         # TODO: to be moved into backend (and attributes?)
 
         prob_zero_events = self.formulae.trivia.poissonian_avoidance_function
 
-        def backend_deposition_removal_all_or_nothing(
+        def backend_sedimentation_removal_all_or_nothing(
             *, relative_fall_velocity, length_scale, timestep, u01
         ):
             pass
 
-        def backend_deposition_removal_deterministic(
+        def backend_sedimentation_removal_deterministic(
             *, relative_fall_velocity, multiplicity, length_scale, timestep
         ):
             for i, velocity in enumerate(relative_fall_velocity):
@@ -586,14 +586,14 @@ def backend_deposition_removal_deterministic(
                 multiplicity[i] *= survive_prob
 
         if all_or_nothing:
-            backend_deposition_removal_all_or_nothing(
+            backend_sedimentation_removal_all_or_nothing(
                 relative_fall_velocity=self.attributes["relative fall velocity"].data,
                 length_scale=length_scale,
                 timestep=self.dt,
                 u01=None,
             )
         else:
-            backend_deposition_removal_deterministic(
+            backend_sedimentation_removal_deterministic(
                 relative_fall_velocity=self.attributes["relative fall velocity"].data,
                 multiplicity=self.attributes["multiplicity"].data,
                 length_scale=length_scale,
diff --git a/tests/unit_tests/dynamics/test_deposition_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py
similarity index 83%
rename from tests/unit_tests/dynamics/test_deposition_removal.py
rename to tests/unit_tests/dynamics/test_sedimentation_removal.py
index 58d127c027..6bd32083f8 100644
--- a/tests/unit_tests/dynamics/test_deposition_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py
@@ -1,4 +1,4 @@
-from PySDM.dynamics import DepositionRemoval
+from PySDM.dynamics import SedimentationRemoval
 from PySDM.physics import si
 from PySDM.environments import Box
 from PySDM.builder import Builder
@@ -9,10 +9,10 @@
 import numpy as np
 
 
-class TestDepositionRemoval:
+class TestSedimentationRemoval:
     @staticmethod
     @pytest.mark.parametrize("all_or_nothing", (True, False))
-    def test_convergence_wrt_dt(all_or_nothing, plot=True):
+    def test_convergence_wrt_dt(all_or_nothing, plot=False):
         # arrange
         dt = 1 * si.s
         dv = 666 * si.m**3
@@ -26,7 +26,7 @@ def test_convergence_wrt_dt(all_or_nothing, plot=True):
             environment=Box(dv=dv, dt=dt),
             backend=backend_instance,
         )
-        builder.add_dynamic(DepositionRemoval(all_or_nothing=all_or_nothing))
+        builder.add_dynamic(SedimentationRemoval(all_or_nothing=all_or_nothing))
         particulator = builder.build(
             attributes={
                 "multiplicity": np.asarray(multiplicity),
@@ -47,10 +47,12 @@ def test_convergence_wrt_dt(all_or_nothing, plot=True):
         pyplot.title(f"{all_or_nothing=}")
         pyplot.xlabel("time [s]")
         pyplot.ylabel("particle concentration [m$^{-3}$]")
-        pyplot.plot(output["time"], output["particle concentration"])
+        pyplot.semilogy(output["time"], output["particle concentration"])
+
         if plot:
             pyplot.show()
         else:
             pyplot.clf()
 
         # assert
+        # TODO!

From f05ae89753586b8e1c9ad73cde262d6375b7f9da Mon Sep 17 00:00:00 2001
From: Sylwester Arabas <sylwester.arabas@agh.edu.pl>
Date: Sat, 8 Nov 2025 15:34:46 +0100
Subject: [PATCH 03/18] extended test cases

---
 .../dynamics/test_sedimentation_removal.py    | 75 ++++++++++++-------
 1 file changed, 47 insertions(+), 28 deletions(-)

diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py
index 6bd32083f8..a915f6bf39 100644
--- a/tests/unit_tests/dynamics/test_sedimentation_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py
@@ -11,43 +11,62 @@
 
 class TestSedimentationRemoval:
     @staticmethod
-    @pytest.mark.parametrize("all_or_nothing", (True, False))
-    def test_convergence_wrt_dt(all_or_nothing, plot=False):
+    @pytest.mark.parametrize("all_or_nothing", (False,))
+    def test_convergence_wrt_dt(all_or_nothing, plot=True):
         # arrange
-        dt = 1 * si.s
-        dv = 666 * si.m**3
-        n_steps = 100
-        multiplicity = [1, 10, 100, 1000]
-        water_mass = [1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug]
+        dts = 5 * si.s, 0.5 * si.s
+        dvs = 1e2 * si.m**3, 1e3 * si.m**3, 1e4 * si.m**3
+        t_max = 500 * si.s
+        multiplicities = 1e5, 1e6, 1e7, 1e8
+        water_masses = 1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug
         backend_instance = CPU()
 
-        builder = Builder(
-            n_sd=len(multiplicity),
-            environment=Box(dv=dv, dt=dt),
-            backend=backend_instance,
-        )
-        builder.add_dynamic(SedimentationRemoval(all_or_nothing=all_or_nothing))
-        particulator = builder.build(
-            attributes={
-                "multiplicity": np.asarray(multiplicity),
-                "signed water mass": np.asarray(water_mass),
-            },
-            products=(ParticleConcentration(), SuperDropletCountPerGridbox(), Time()),
-        )
-
         # act
-        output = {name: [] for name in particulator.products}
-        for step in range(n_steps):
-            if step != 0:
-                particulator.run(steps=1)
-            for name, product in particulator.products.items():
-                output[name].append(product.get() + 0)
+        output = {}
+        for dt in dts:
+            for dv in dvs:
+                builder = Builder(
+                    n_sd=len(multiplicities),
+                    environment=Box(dv=dv, dt=dt),
+                    backend=backend_instance,
+                )
+                builder.add_dynamic(SedimentationRemoval(all_or_nothing=all_or_nothing))
+                particulator = builder.build(
+                    attributes={
+                        "multiplicity": np.asarray(multiplicities),
+                        "signed water mass": np.asarray(water_masses),
+                    },
+                    products=(
+                        ParticleConcentration(),
+                        SuperDropletCountPerGridbox(),
+                        Time(),
+                    ),
+                )
+                key = f"{dt=} {dv=}"
+                output[key] = {name: [] for name in particulator.products}
+                while particulator.n_steps * dt <= t_max:
+                    if len(output[key]["time"]) != 0:
+                        particulator.run(steps=1)
+                    for name, product in particulator.products.items():
+                        output[key][name].append(product.get() + 0)
 
         # plot
         pyplot.title(f"{all_or_nothing=}")
         pyplot.xlabel("time [s]")
         pyplot.ylabel("particle concentration [m$^{-3}$]")
-        pyplot.semilogy(output["time"], output["particle concentration"])
+        for dt in dts:
+            for dv in dvs:
+                key = f"{dt=} {dv=}"
+                pyplot.plot(
+                    output[key]["time"],
+                    output[key]["particle concentration"],
+                    label=key,
+                    linewidth=4 + 3 * np.log10(dt),
+                )
+        pyplot.gca().set_yscale("log")
+        pyplot.gca().set_xlim(left=0, right=t_max)
+        pyplot.legend()
+        pyplot.grid()
 
         if plot:
             pyplot.show()

From b253d093a3664cfc312bd8b954a69a858e5103e2 Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Wed, 20 May 2026 12:32:39 +0200
Subject: [PATCH 04/18] renaming functions and booleans

---
 PySDM/dynamics/sedimentation_removal.py            |  6 +++---
 PySDM/particulator.py                              | 14 +++++++-------
 .../dynamics/test_sedimentation_removal.py         | 10 ++++++----
 3 files changed, 16 insertions(+), 14 deletions(-)

diff --git a/PySDM/dynamics/sedimentation_removal.py b/PySDM/dynamics/sedimentation_removal.py
index bdb3414f7e..3b1e50cfb2 100644
--- a/PySDM/dynamics/sedimentation_removal.py
+++ b/PySDM/dynamics/sedimentation_removal.py
@@ -6,9 +6,9 @@
 
 @register_dynamic()
 class SedimentationRemoval:
-    def __init__(self, *, all_or_nothing: bool):
+    def __init__(self, *, stochastic_deposition_removal: bool):
         """stochastic ("all or nothing") or deterministic (multiplicity altering) removal"""
-        self.all_or_nothing = all_or_nothing
+        self.stochastic_deposition_removal = stochastic_deposition_removal
         self.length_scale = None
 
     def register(self, builder):
@@ -21,6 +21,6 @@ def __call__(self):
         """see, e.g., the naive scheme in Algorithm 1 in
         [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)"""
         self.particulator.sedimentation_removal(
-            all_or_nothing=self.all_or_nothing,
+            stochastic_deposition_removal=self.stochastic_deposition_removal,
             length_scale=self.length_scale,
         )
diff --git a/PySDM/particulator.py b/PySDM/particulator.py
index c43e0ec740..f1b6423b69 100644
--- a/PySDM/particulator.py
+++ b/PySDM/particulator.py
@@ -593,17 +593,17 @@ def thaw_instantaneous(self):
             temperature=self.environment["T"],
         )
 
-    def sedimentation_removal(self, *, all_or_nothing, length_scale):
+    def sedimentation_removal(self, *, stochastic_deposition_removal, length_scale):
         # TODO: to be moved into backend (and attributes?)
 
         prob_zero_events = self.formulae.trivia.poissonian_avoidance_function
 
-        def backend_sedimentation_removal_all_or_nothing(
-            *, relative_fall_velocity, length_scale, timestep, u01
+        def backend_sedimentation_removal_deterministic(
+            *, relative_fall_velocity, multiplicity, length_scale, timestep
         ):
             pass
 
-        def backend_sedimentation_removal_deterministic(
+        def backend_sedimentation_removal_stochastic(
             *, relative_fall_velocity, multiplicity, length_scale, timestep
         ):
             for i, velocity in enumerate(relative_fall_velocity):
@@ -612,12 +612,12 @@ def backend_sedimentation_removal_deterministic(
                 assert 0 <= survive_prob <= 1
                 multiplicity[i] *= survive_prob
 
-        if all_or_nothing:
-            backend_sedimentation_removal_all_or_nothing(
+        if stochastic_deposition_removal:
+            backend_sedimentation_removal_stochastic(
                 relative_fall_velocity=self.attributes["relative fall velocity"].data,
+                multiplicity=self.attributes["multiplicity"].data,
                 length_scale=length_scale,
                 timestep=self.dt,
-                u01=None,
             )
         else:
             backend_sedimentation_removal_deterministic(
diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py
index a915f6bf39..7ced62d1e6 100644
--- a/tests/unit_tests/dynamics/test_sedimentation_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py
@@ -11,8 +11,8 @@
 
 class TestSedimentationRemoval:
     @staticmethod
-    @pytest.mark.parametrize("all_or_nothing", (False,))
-    def test_convergence_wrt_dt(all_or_nothing, plot=True):
+    @pytest.mark.parametrize("stochastic", (True,))
+    def test_convergence_wrt_dt(stochastic, plot=True):
         # arrange
         dts = 5 * si.s, 0.5 * si.s
         dvs = 1e2 * si.m**3, 1e3 * si.m**3, 1e4 * si.m**3
@@ -30,7 +30,9 @@ def test_convergence_wrt_dt(all_or_nothing, plot=True):
                     environment=Box(dv=dv, dt=dt),
                     backend=backend_instance,
                 )
-                builder.add_dynamic(SedimentationRemoval(all_or_nothing=all_or_nothing))
+                builder.add_dynamic(
+                    SedimentationRemoval(stochastic_deposition_removal=stochastic)
+                )
                 particulator = builder.build(
                     attributes={
                         "multiplicity": np.asarray(multiplicities),
@@ -51,7 +53,7 @@ def test_convergence_wrt_dt(all_or_nothing, plot=True):
                         output[key][name].append(product.get() + 0)
 
         # plot
-        pyplot.title(f"{all_or_nothing=}")
+        pyplot.title(f"{stochastic=}")
         pyplot.xlabel("time [s]")
         pyplot.ylabel("particle concentration [m$^{-3}$]")
         for dt in dts:

From 15530fbbf679ca88bb9ddfeb080a9064360e8c0f Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Wed, 20 May 2026 14:03:44 +0200
Subject: [PATCH 05/18] moved sedimentation_removal to backend

---
 PySDM/backends/impl_numba/methods/__init__.py |  1 +
 PySDM/backends/numba.py                       |  1 +
 PySDM/particulator.py                         | 20 ++-----------------
 3 files changed, 4 insertions(+), 18 deletions(-)

diff --git a/PySDM/backends/impl_numba/methods/__init__.py b/PySDM/backends/impl_numba/methods/__init__.py
index 4c613462a2..808be56292 100644
--- a/PySDM/backends/impl_numba/methods/__init__.py
+++ b/PySDM/backends/impl_numba/methods/__init__.py
@@ -14,3 +14,4 @@
 from .terminal_velocity_methods import TerminalVelocityMethods
 from .seeding_methods import SeedingMethods
 from .deposition_methods import DepositionMethods
+from .sedimentation_removal_methods import SedimentationRemovalMethods
diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py
index 21f82478ed..32d2992625 100644
--- a/PySDM/backends/numba.py
+++ b/PySDM/backends/numba.py
@@ -32,6 +32,7 @@ class Numba(  # pylint: disable=too-many-ancestors,duplicate-code
     methods.IsotopeMethods,
     methods.SeedingMethods,
     methods.DepositionMethods,
+    methods.SedimentationRemovalMethods,
 ):
     Storage = ImportedStorage
     Random = ImportedRandom
diff --git a/PySDM/particulator.py b/PySDM/particulator.py
index f1b6423b69..246789214e 100644
--- a/PySDM/particulator.py
+++ b/PySDM/particulator.py
@@ -596,31 +596,15 @@ def thaw_instantaneous(self):
     def sedimentation_removal(self, *, stochastic_deposition_removal, length_scale):
         # TODO: to be moved into backend (and attributes?)
 
-        prob_zero_events = self.formulae.trivia.poissonian_avoidance_function
-
-        def backend_sedimentation_removal_deterministic(
-            *, relative_fall_velocity, multiplicity, length_scale, timestep
-        ):
-            pass
-
-        def backend_sedimentation_removal_stochastic(
-            *, relative_fall_velocity, multiplicity, length_scale, timestep
-        ):
-            for i, velocity in enumerate(relative_fall_velocity):
-                removal_rate = velocity / length_scale
-                survive_prob = prob_zero_events(r=removal_rate, dt=timestep)
-                assert 0 <= survive_prob <= 1
-                multiplicity[i] *= survive_prob
-
         if stochastic_deposition_removal:
-            backend_sedimentation_removal_stochastic(
+            self.backend.sedimentation_removal_stochastic(
                 relative_fall_velocity=self.attributes["relative fall velocity"].data,
                 multiplicity=self.attributes["multiplicity"].data,
                 length_scale=length_scale,
                 timestep=self.dt,
             )
         else:
-            backend_sedimentation_removal_deterministic(
+            self.backend.sedimentation_removal_deterministic(
                 relative_fall_velocity=self.attributes["relative fall velocity"].data,
                 multiplicity=self.attributes["multiplicity"].data,
                 length_scale=length_scale,

From 073fc48f4bff403706eb511b8092b8ce6f214729 Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Wed, 20 May 2026 14:25:28 +0200
Subject: [PATCH 06/18] Added jit commands to new method. Included ice in unit
 testing

---
 .../methods/sedimentation_removal_methods.py  | 51 +++++++++++++++++++
 .../dynamics/test_sedimentation_removal.py    |  4 +-
 2 files changed, 53 insertions(+), 2 deletions(-)
 create mode 100644 PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py

diff --git a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
new file mode 100644
index 0000000000..e037da0d0b
--- /dev/null
+++ b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
@@ -0,0 +1,51 @@
+"""
+CPU implementation of backend methods for removal due to sedimentation
+in a 0D environment
+"""
+
+from functools import cached_property
+
+import numba
+
+from PySDM.backends.impl_common.backend_methods import BackendMethods
+
+
+class SedimentationRemovalMethods(BackendMethods):
+
+    @cached_property
+    def _sedimentation_removal_deterministic_body(self):
+
+        @numba.njit(**self.default_jit_flags)
+        def body(relative_fall_velocity, multiplicity, length_scale, timestep):
+            pass
+
+        return body
+
+    @cached_property
+    def _sedimentation_removal_stochastic_body(self):
+
+        prob_zero_events = self.formulae.trivia.poissonian_avoidance_function
+
+        @numba.njit(**self.default_jit_flags)
+        def body(relative_fall_velocity, multiplicity, length_scale, timestep):
+            for i, velocity in enumerate(relative_fall_velocity):
+                removal_rate = velocity / length_scale
+                survive_prob = prob_zero_events(r=removal_rate, dt=timestep)
+                assert 0 <= survive_prob <= 1
+                multiplicity[i] *= survive_prob
+
+        return body
+
+    def sedimentation_removal_deterministic(
+        self, *, relative_fall_velocity, multiplicity, length_scale, timestep
+    ):
+        self._sedimentation_removal_deterministic_body(
+            relative_fall_velocity, multiplicity, length_scale, timestep
+        )
+
+    def sedimentation_removal_stochastic(
+        self, *, relative_fall_velocity, multiplicity, length_scale, timestep
+    ):
+        self._sedimentation_removal_stochastic_body(
+            relative_fall_velocity, multiplicity, length_scale, timestep
+        )
diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py
index 7ced62d1e6..0f1277af77 100644
--- a/tests/unit_tests/dynamics/test_sedimentation_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py
@@ -17,8 +17,8 @@ def test_convergence_wrt_dt(stochastic, plot=True):
         dts = 5 * si.s, 0.5 * si.s
         dvs = 1e2 * si.m**3, 1e3 * si.m**3, 1e4 * si.m**3
         t_max = 500 * si.s
-        multiplicities = 1e5, 1e6, 1e7, 1e8
-        water_masses = 1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug
+        multiplicities = 1e5, 1e6, 1e7, 1e8, 1e5
+        water_masses = 1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug, -1 * si.ug
         backend_instance = CPU()
 
         # act

From a9815371cccc39a6940365ce300501cc883ed3ab Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Wed, 20 May 2026 15:29:30 +0200
Subject: [PATCH 07/18] Added determernistic sedimentation removal

---
 .../impl_numba/methods/sedimentation_removal_methods.py    | 3 ++-
 PySDM/dynamics/sedimentation_removal.py                    | 7 ++++---
 PySDM/particulator.py                                      | 2 --
 tests/unit_tests/dynamics/test_sedimentation_removal.py    | 2 +-
 4 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
index e037da0d0b..22c79b2f18 100644
--- a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
+++ b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
@@ -17,7 +17,8 @@ def _sedimentation_removal_deterministic_body(self):
 
         @numba.njit(**self.default_jit_flags)
         def body(relative_fall_velocity, multiplicity, length_scale, timestep):
-            pass
+            for i, velocity in enumerate(relative_fall_velocity):
+                multiplicity[i] -= multiplicity[i] * velocity * timestep / length_scale
 
         return body
 
diff --git a/PySDM/dynamics/sedimentation_removal.py b/PySDM/dynamics/sedimentation_removal.py
index 3b1e50cfb2..6ad4e95bed 100644
--- a/PySDM/dynamics/sedimentation_removal.py
+++ b/PySDM/dynamics/sedimentation_removal.py
@@ -1,13 +1,14 @@
 """deposition removal logic for zero-dimensional environments"""
 
+from typing import Optional
 from PySDM.dynamics.impl import register_dynamic
 import numpy as np
 
 
 @register_dynamic()
 class SedimentationRemoval:
-    def __init__(self, *, stochastic_deposition_removal: bool):
-        """stochastic ("all or nothing") or deterministic (multiplicity altering) removal"""
+    def __init__(self, *, stochastic_deposition_removal: Optional[bool] = True):
+        """stochastic or deterministic removal"""
         self.stochastic_deposition_removal = stochastic_deposition_removal
         self.length_scale = None
 
@@ -18,7 +19,7 @@ def register(self, builder):
         self.length_scale = np.cbrt(self.particulator.environment.mesh.dv)
 
     def __call__(self):
-        """see, e.g., the naive scheme in Algorithm 1 in
+        """for stochastic removal, see, e.g., the naive scheme in Algorithm 1 in
         [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)"""
         self.particulator.sedimentation_removal(
             stochastic_deposition_removal=self.stochastic_deposition_removal,
diff --git a/PySDM/particulator.py b/PySDM/particulator.py
index 246789214e..95fb9e4702 100644
--- a/PySDM/particulator.py
+++ b/PySDM/particulator.py
@@ -594,8 +594,6 @@ def thaw_instantaneous(self):
         )
 
     def sedimentation_removal(self, *, stochastic_deposition_removal, length_scale):
-        # TODO: to be moved into backend (and attributes?)
-
         if stochastic_deposition_removal:
             self.backend.sedimentation_removal_stochastic(
                 relative_fall_velocity=self.attributes["relative fall velocity"].data,
diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py
index 0f1277af77..13ae6a3198 100644
--- a/tests/unit_tests/dynamics/test_sedimentation_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py
@@ -11,7 +11,7 @@
 
 class TestSedimentationRemoval:
     @staticmethod
-    @pytest.mark.parametrize("stochastic", (True,))
+    @pytest.mark.parametrize("stochastic", (True, False))
     def test_convergence_wrt_dt(stochastic, plot=True):
         # arrange
         dts = 5 * si.s, 0.5 * si.s

From edf5766758a4f487d430dfbfc5395e4f78b01cc6 Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Wed, 20 May 2026 17:46:48 +0200
Subject: [PATCH 08/18] Added asserts to unit test

---
 .../dynamics/test_sedimentation_removal.py    | 30 ++++++++++++++++---
 1 file changed, 26 insertions(+), 4 deletions(-)

diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py
index 13ae6a3198..0d35af02c0 100644
--- a/tests/unit_tests/dynamics/test_sedimentation_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py
@@ -12,11 +12,11 @@
 class TestSedimentationRemoval:
     @staticmethod
     @pytest.mark.parametrize("stochastic", (True, False))
-    def test_convergence_wrt_dt(stochastic, plot=True):
+    def test_convergence_wrt_dt(stochastic, plot=False):
         # arrange
-        dts = 5 * si.s, 0.5 * si.s
+        dts = 0.5 * si.s, 5 * si.s
         dvs = 1e2 * si.m**3, 1e3 * si.m**3, 1e4 * si.m**3
-        t_max = 500 * si.s
+        t_max = 600 * si.s
         multiplicities = 1e5, 1e6, 1e7, 1e8, 1e5
         water_masses = 1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug, -1 * si.ug
         backend_instance = CPU()
@@ -76,4 +76,26 @@ def test_convergence_wrt_dt(stochastic, plot=True):
             pyplot.clf()
 
         # assert
-        # TODO!
+        for dt in dts:
+            for dv in dvs:
+                key = f"{dt=} {dv=}"
+                particle_concentration = np.asarray(
+                    output[key]["particle concentration"]
+                )
+                assert particle_concentration[-1] == 0.0
+
+        for dv in dvs:
+            time_compare = 0
+            for dt in dts:
+                key = f"{dt=} {dv=}"
+                time_end = np.asarray(output[key]["time"])[-1]
+                assert time_end >= time_compare
+                time_compare = time_end
+
+        for dt in dts:
+            time_compare = 0
+            for dv in dvs:
+                key = f"{dt=} {dv=}"
+                time_end = np.asarray(output[key]["time"])[-1]
+                assert time_end >= time_compare
+                time_compare = time_end

From 8584ef53e111766fcf219081669ba4caf0b744f6 Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Thu, 21 May 2026 12:52:21 +0200
Subject: [PATCH 09/18] pylint and pdoc and unit-tests fix

---
 PySDM/dynamics/__init__.py                    |  2 +-
 PySDM/dynamics/sedimentation_removal.py       | 11 ++++---
 PySDM/particulator.py                         |  4 +--
 docs/bibliography.json                        |  2 +-
 .../dynamics/test_sedimentation_removal.py    | 31 +++++++++++--------
 5 files changed, 28 insertions(+), 22 deletions(-)

diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py
index 45725a2dac..90b107870b 100644
--- a/PySDM/dynamics/__init__.py
+++ b/PySDM/dynamics/__init__.py
@@ -17,4 +17,4 @@
 from PySDM.dynamics.relaxed_velocity import RelaxedVelocity
 from PySDM.dynamics.seeding import Seeding
 from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce
-from PySDM.dynamics.sedimentation_removal import SedimentationRemoval
+from PySDM.dynamics.sedimentation_removal import ParcelEnvironmentSedimentationRemoval
diff --git a/PySDM/dynamics/sedimentation_removal.py b/PySDM/dynamics/sedimentation_removal.py
index 6ad4e95bed..9f39e2a484 100644
--- a/PySDM/dynamics/sedimentation_removal.py
+++ b/PySDM/dynamics/sedimentation_removal.py
@@ -1,16 +1,17 @@
 """deposition removal logic for zero-dimensional environments"""
 
+import numpy as np
 from typing import Optional
 from PySDM.dynamics.impl import register_dynamic
-import numpy as np
 
 
 @register_dynamic()
-class SedimentationRemoval:
-    def __init__(self, *, stochastic_deposition_removal: Optional[bool] = True):
+class ParcelEnvironmentSedimentationRemoval:
+    def __init__(self, *, stochastic_sedimentation_removal: Optional[bool] = True):
         """stochastic or deterministic removal"""
-        self.stochastic_deposition_removal = stochastic_deposition_removal
+        self.stochastic_sedimentation_removal = stochastic_sedimentation_removal
         self.length_scale = None
+        self.particulator = None
 
     def register(self, builder):
         builder.request_attribute("relative fall velocity")
@@ -22,6 +23,6 @@ def __call__(self):
         """for stochastic removal, see, e.g., the naive scheme in Algorithm 1 in
         [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)"""
         self.particulator.sedimentation_removal(
-            stochastic_deposition_removal=self.stochastic_deposition_removal,
+            stochastic_sedimentation_removal=self.stochastic_sedimentation_removal,
             length_scale=self.length_scale,
         )
diff --git a/PySDM/particulator.py b/PySDM/particulator.py
index 95fb9e4702..3257033d52 100644
--- a/PySDM/particulator.py
+++ b/PySDM/particulator.py
@@ -593,8 +593,8 @@ def thaw_instantaneous(self):
             temperature=self.environment["T"],
         )
 
-    def sedimentation_removal(self, *, stochastic_deposition_removal, length_scale):
-        if stochastic_deposition_removal:
+    def sedimentation_removal(self, *, stochastic_sedimentation_removal, length_scale):
+        if stochastic_sedimentation_removal:
             self.backend.sedimentation_removal_stochastic(
                 relative_fall_velocity=self.attributes["relative fall velocity"].data,
                 multiplicity=self.attributes["multiplicity"].data,
diff --git a/docs/bibliography.json b/docs/bibliography.json
index 6f28515ab0..94280a63cc 100644
--- a/docs/bibliography.json
+++ b/docs/bibliography.json
@@ -1041,7 +1041,7 @@
   },
   "https://doi.org/10.1016/j.jcp.2016.06.029": {
     "usages": [
-      "PySDM/dynamics/deposition_removal.py"
+      "PySDM/dynamics/sedimentation_removal.py"
     ],
     "label": "Curtis et al. 2016 (J. Comp. Phys. 322)",
     "title": "Accelerated simulation of stochastic particle removal\nprocesses in particle-resolved aerosol models"
diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py
index 0d35af02c0..e246726bb9 100644
--- a/tests/unit_tests/dynamics/test_sedimentation_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py
@@ -1,25 +1,28 @@
-from PySDM.dynamics import SedimentationRemoval
-from PySDM.physics import si
-from PySDM.environments import Box
-from PySDM.builder import Builder
-from PySDM.backends import CPU
-from PySDM.products import ParticleConcentration, SuperDropletCountPerGridbox, Time
+# pylint: disable=missing-module-docstring
 from matplotlib import pyplot
 import pytest
 import numpy as np
+from PySDM.dynamics import ParcelEnvironmentSedimentationRemoval
+from PySDM.physics import si
+from PySDM.environments import Box
+from PySDM import Builder, Formulae
+from PySDM.backends import ThrustRTC
+from PySDM.products import ParticleConcentration, SuperDropletCountPerGridbox, Time
 
 
-class TestSedimentationRemoval:
+class TestSedimentationRemoval:  # pylint: disable=too-few-public-methods,too-many-branches,too-many-locals
     @staticmethod
     @pytest.mark.parametrize("stochastic", (True, False))
-    def test_convergence_wrt_dt(stochastic, plot=False):
+    def test_convergence_wrt_dt(backend_class, stochastic, plot=False):
         # arrange
         dts = 0.5 * si.s, 5 * si.s
         dvs = 1e2 * si.m**3, 1e3 * si.m**3, 1e4 * si.m**3
         t_max = 600 * si.s
-        multiplicities = 1e5, 1e6, 1e7, 1e8, 1e5
-        water_masses = 1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug, -1 * si.ug
-        backend_instance = CPU()
+        multiplicities = 1e5, 1e6, 1e7, 1e8
+        water_masses = 1 * si.ug, 2 * si.ug, 3 * si.ug, 4 * si.ug
+
+        if backend_class is ThrustRTC:
+            pytest.skip("TODO #1871")
 
         # act
         output = {}
@@ -28,10 +31,12 @@ def test_convergence_wrt_dt(stochastic, plot=False):
                 builder = Builder(
                     n_sd=len(multiplicities),
                     environment=Box(dv=dv, dt=dt),
-                    backend=backend_instance,
+                    backend=backend_class(),
                 )
                 builder.add_dynamic(
-                    SedimentationRemoval(stochastic_deposition_removal=stochastic)
+                    ParcelEnvironmentSedimentationRemoval(
+                        stochastic_sedimentation_removal=stochastic
+                    )
                 )
                 particulator = builder.build(
                     attributes={

From 442dc5babbd503f7a862593dd7565891ba701c6b Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Thu, 21 May 2026 13:52:55 +0200
Subject: [PATCH 10/18] addressing pylint

---
 PySDM/dynamics/sedimentation_removal.py                 | 2 +-
 tests/unit_tests/dynamics/test_sedimentation_removal.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/PySDM/dynamics/sedimentation_removal.py b/PySDM/dynamics/sedimentation_removal.py
index 9f39e2a484..5ad48eaa9b 100644
--- a/PySDM/dynamics/sedimentation_removal.py
+++ b/PySDM/dynamics/sedimentation_removal.py
@@ -1,7 +1,7 @@
 """deposition removal logic for zero-dimensional environments"""
 
-import numpy as np
 from typing import Optional
+import numpy as np
 from PySDM.dynamics.impl import register_dynamic
 
 
diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal.py
index e246726bb9..ba9f7fdb90 100644
--- a/tests/unit_tests/dynamics/test_sedimentation_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal.py
@@ -5,7 +5,7 @@
 from PySDM.dynamics import ParcelEnvironmentSedimentationRemoval
 from PySDM.physics import si
 from PySDM.environments import Box
-from PySDM import Builder, Formulae
+from PySDM import Builder
 from PySDM.backends import ThrustRTC
 from PySDM.products import ParticleConcentration, SuperDropletCountPerGridbox, Time
 

From 224e3601de9da49379f6e84893c1308336aedf04 Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Thu, 21 May 2026 15:12:49 +0200
Subject: [PATCH 11/18] numba fix?

---
 .../methods/sedimentation_removal_methods.py      | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
index 22c79b2f18..0e88a872be 100644
--- a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
+++ b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
@@ -17,8 +17,14 @@ def _sedimentation_removal_deterministic_body(self):
 
         @numba.njit(**self.default_jit_flags)
         def body(relative_fall_velocity, multiplicity, length_scale, timestep):
-            for i, velocity in enumerate(relative_fall_velocity):
-                multiplicity[i] -= multiplicity[i] * velocity * timestep / length_scale
+            n_sd = len(relative_fall_velocity)
+            for i in numba.prange(n_sd):
+                multiplicity[i] -= (
+                    multiplicity[i]
+                    * relative_fall_velocity[i]
+                    * timestep
+                    / length_scale
+                )
 
         return body
 
@@ -29,8 +35,9 @@ def _sedimentation_removal_stochastic_body(self):
 
         @numba.njit(**self.default_jit_flags)
         def body(relative_fall_velocity, multiplicity, length_scale, timestep):
-            for i, velocity in enumerate(relative_fall_velocity):
-                removal_rate = velocity / length_scale
+            n_sd = len(relative_fall_velocity)
+            for i in numba.prange(n_sd):
+                removal_rate = relative_fall_velocity[i] / length_scale
                 survive_prob = prob_zero_events(r=removal_rate, dt=timestep)
                 assert 0 <= survive_prob <= 1
                 multiplicity[i] *= survive_prob

From 889d85eba05cbe5eaf730f382d18efd88e61806e Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Thu, 21 May 2026 15:45:44 +0200
Subject: [PATCH 12/18] added pylint comment

---
 .../impl_numba/methods/sedimentation_removal_methods.py       | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
index 0e88a872be..a7485d14bf 100644
--- a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
+++ b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
@@ -18,7 +18,7 @@ def _sedimentation_removal_deterministic_body(self):
         @numba.njit(**self.default_jit_flags)
         def body(relative_fall_velocity, multiplicity, length_scale, timestep):
             n_sd = len(relative_fall_velocity)
-            for i in numba.prange(n_sd):
+            for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
                 multiplicity[i] -= (
                     multiplicity[i]
                     * relative_fall_velocity[i]
@@ -36,7 +36,7 @@ def _sedimentation_removal_stochastic_body(self):
         @numba.njit(**self.default_jit_flags)
         def body(relative_fall_velocity, multiplicity, length_scale, timestep):
             n_sd = len(relative_fall_velocity)
-            for i in numba.prange(n_sd):
+            for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
                 removal_rate = relative_fall_velocity[i] / length_scale
                 survive_prob = prob_zero_events(r=removal_rate, dt=timestep)
                 assert 0 <= survive_prob <= 1

From 31072a2517042591ed2c5e3f30613ecac78e421e Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Wed, 27 May 2026 11:54:55 +0200
Subject: [PATCH 13/18] Removed assertion from loop

---
 .../backends/impl_numba/methods/sedimentation_removal_methods.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
index a7485d14bf..9955fa48d4 100644
--- a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
+++ b/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
@@ -39,7 +39,6 @@ def body(relative_fall_velocity, multiplicity, length_scale, timestep):
             for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
                 removal_rate = relative_fall_velocity[i] / length_scale
                 survive_prob = prob_zero_events(r=removal_rate, dt=timestep)
-                assert 0 <= survive_prob <= 1
                 multiplicity[i] *= survive_prob
 
         return body

From f37742d1e2fb9b4f17ecb1b7beae700dbb27f715 Mon Sep 17 00:00:00 2001
From: Sylwester Arabas <sylwester.arabas@agh.edu.pl>
Date: Sat, 30 May 2026 16:30:12 +0200
Subject: [PATCH 14/18] remove spurious "\n" from paper title in bib JSON

---
 docs/bibliography.json | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/bibliography.json b/docs/bibliography.json
index e4e6a690a3..3f4ab8eb71 100644
--- a/docs/bibliography.json
+++ b/docs/bibliography.json
@@ -1053,6 +1053,6 @@
       "PySDM/dynamics/sedimentation_removal.py"
     ],
     "label": "Curtis et al. 2016 (J. Comp. Phys. 322)",
-    "title": "Accelerated simulation of stochastic particle removal\nprocesses in particle-resolved aerosol models"
+    "title": "Accelerated simulation of stochastic particle removal processes in particle-resolved aerosol models"
   }
 }

From 2d3a243b168c0e98b3d9e079c35acc041e2e0480 Mon Sep 17 00:00:00 2001
From: Sylwester Arabas <sylwester.arabas@agh.edu.pl>
Date: Sat, 30 May 2026 16:47:13 +0200
Subject: [PATCH 15/18] fix compatibility with newer PySDM API (no
 add_dynamic); rename the new dynamic to SedimentationRemoval0D to hint it can
 also be used with Box environment

---
 PySDM/backends/impl_numba/methods/__init__.py      |  2 +-
 ...hods.py => sedimentation_removal_0d_methods.py} |  2 +-
 PySDM/backends/numba.py                            |  2 +-
 PySDM/dynamics/__init__.py                         |  2 +-
 ...tion_removal.py => sedimentation_removal_0d.py} |  2 +-
 ...removal.py => test_sedimentation_removal_0d.py} | 14 +++++++-------
 6 files changed, 12 insertions(+), 12 deletions(-)
 rename PySDM/backends/impl_numba/methods/{sedimentation_removal_methods.py => sedimentation_removal_0d_methods.py} (97%)
 rename PySDM/dynamics/{sedimentation_removal.py => sedimentation_removal_0d.py} (96%)
 rename tests/unit_tests/dynamics/{test_sedimentation_removal.py => test_sedimentation_removal_0d.py} (90%)

diff --git a/PySDM/backends/impl_numba/methods/__init__.py b/PySDM/backends/impl_numba/methods/__init__.py
index 808be56292..192cc0c6c7 100644
--- a/PySDM/backends/impl_numba/methods/__init__.py
+++ b/PySDM/backends/impl_numba/methods/__init__.py
@@ -14,4 +14,4 @@
 from .terminal_velocity_methods import TerminalVelocityMethods
 from .seeding_methods import SeedingMethods
 from .deposition_methods import DepositionMethods
-from .sedimentation_removal_methods import SedimentationRemovalMethods
+from .sedimentation_removal_0d_methods import SedimentationRemoval0DMethods
diff --git a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py b/PySDM/backends/impl_numba/methods/sedimentation_removal_0d_methods.py
similarity index 97%
rename from PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
rename to PySDM/backends/impl_numba/methods/sedimentation_removal_0d_methods.py
index 9955fa48d4..2c090af409 100644
--- a/PySDM/backends/impl_numba/methods/sedimentation_removal_methods.py
+++ b/PySDM/backends/impl_numba/methods/sedimentation_removal_0d_methods.py
@@ -10,7 +10,7 @@
 from PySDM.backends.impl_common.backend_methods import BackendMethods
 
 
-class SedimentationRemovalMethods(BackendMethods):
+class SedimentationRemoval0DMethods(BackendMethods):
 
     @cached_property
     def _sedimentation_removal_deterministic_body(self):
diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py
index 32d2992625..a7826c968b 100644
--- a/PySDM/backends/numba.py
+++ b/PySDM/backends/numba.py
@@ -32,7 +32,7 @@ class Numba(  # pylint: disable=too-many-ancestors,duplicate-code
     methods.IsotopeMethods,
     methods.SeedingMethods,
     methods.DepositionMethods,
-    methods.SedimentationRemovalMethods,
+    methods.SedimentationRemoval0DMethods,
 ):
     Storage = ImportedStorage
     Random = ImportedRandom
diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py
index 90b107870b..11d039fa7a 100644
--- a/PySDM/dynamics/__init__.py
+++ b/PySDM/dynamics/__init__.py
@@ -17,4 +17,4 @@
 from PySDM.dynamics.relaxed_velocity import RelaxedVelocity
 from PySDM.dynamics.seeding import Seeding
 from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce
-from PySDM.dynamics.sedimentation_removal import ParcelEnvironmentSedimentationRemoval
+from PySDM.dynamics.sedimentation_removal_0d import SedimentationRemoval0D
diff --git a/PySDM/dynamics/sedimentation_removal.py b/PySDM/dynamics/sedimentation_removal_0d.py
similarity index 96%
rename from PySDM/dynamics/sedimentation_removal.py
rename to PySDM/dynamics/sedimentation_removal_0d.py
index 5ad48eaa9b..fb79c29e59 100644
--- a/PySDM/dynamics/sedimentation_removal.py
+++ b/PySDM/dynamics/sedimentation_removal_0d.py
@@ -6,7 +6,7 @@
 
 
 @register_dynamic()
-class ParcelEnvironmentSedimentationRemoval:
+class SedimentationRemoval0D:
     def __init__(self, *, stochastic_sedimentation_removal: Optional[bool] = True):
         """stochastic or deterministic removal"""
         self.stochastic_sedimentation_removal = stochastic_sedimentation_removal
diff --git a/tests/unit_tests/dynamics/test_sedimentation_removal.py b/tests/unit_tests/dynamics/test_sedimentation_removal_0d.py
similarity index 90%
rename from tests/unit_tests/dynamics/test_sedimentation_removal.py
rename to tests/unit_tests/dynamics/test_sedimentation_removal_0d.py
index ba9f7fdb90..c8b90e8137 100644
--- a/tests/unit_tests/dynamics/test_sedimentation_removal.py
+++ b/tests/unit_tests/dynamics/test_sedimentation_removal_0d.py
@@ -2,7 +2,7 @@
 from matplotlib import pyplot
 import pytest
 import numpy as np
-from PySDM.dynamics import ParcelEnvironmentSedimentationRemoval
+from PySDM.dynamics import SedimentationRemoval0D
 from PySDM.physics import si
 from PySDM.environments import Box
 from PySDM import Builder
@@ -10,7 +10,7 @@
 from PySDM.products import ParticleConcentration, SuperDropletCountPerGridbox, Time
 
 
-class TestSedimentationRemoval:  # pylint: disable=too-few-public-methods,too-many-branches,too-many-locals
+class TestSedimentationRemoval0D:  # pylint: disable=too-few-public-methods,too-many-branches,too-many-locals
     @staticmethod
     @pytest.mark.parametrize("stochastic", (True, False))
     def test_convergence_wrt_dt(backend_class, stochastic, plot=False):
@@ -32,11 +32,11 @@ def test_convergence_wrt_dt(backend_class, stochastic, plot=False):
                     n_sd=len(multiplicities),
                     environment=Box(dv=dv, dt=dt),
                     backend=backend_class(),
-                )
-                builder.add_dynamic(
-                    ParcelEnvironmentSedimentationRemoval(
-                        stochastic_sedimentation_removal=stochastic
-                    )
+                    dynamics=(
+                        SedimentationRemoval0D(
+                            stochastic_sedimentation_removal=stochastic
+                        ),
+                    ),
                 )
                 particulator = builder.build(
                     attributes={

From aa15e27e92760ceb73a5916a24143ea76097e014 Mon Sep 17 00:00:00 2001
From: Sylwester Arabas <sylwester.arabas@agh.edu.pl>
Date: Sat, 30 May 2026 16:49:17 +0200
Subject: [PATCH 16/18] fix file name in bib

---
 docs/bibliography.json | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/bibliography.json b/docs/bibliography.json
index 3f4ab8eb71..cf8261e379 100644
--- a/docs/bibliography.json
+++ b/docs/bibliography.json
@@ -1050,7 +1050,7 @@
   },
   "https://doi.org/10.1016/j.jcp.2016.06.029": {
     "usages": [
-      "PySDM/dynamics/sedimentation_removal.py"
+      "PySDM/dynamics/sedimentation_removal_0d.py"
     ],
     "label": "Curtis et al. 2016 (J. Comp. Phys. 322)",
     "title": "Accelerated simulation of stochastic particle removal processes in particle-resolved aerosol models"

From 978d341e4c389af26cbfcf1207635b07e2ffa672 Mon Sep 17 00:00:00 2001
From: Tim Luettmer <tluettm@uni-mainz.de>
Date: Mon, 1 Jun 2026 11:50:44 +0200
Subject: [PATCH 17/18] Moved length_scale calculation from register to call

---
 PySDM/dynamics/sedimentation_removal_0d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/PySDM/dynamics/sedimentation_removal_0d.py b/PySDM/dynamics/sedimentation_removal_0d.py
index fb79c29e59..67dfbf8527 100644
--- a/PySDM/dynamics/sedimentation_removal_0d.py
+++ b/PySDM/dynamics/sedimentation_removal_0d.py
@@ -17,11 +17,11 @@ def register(self, builder):
         builder.request_attribute("relative fall velocity")
         assert builder.particulator.environment.mesh.n_dims == 0
         self.particulator = builder.particulator
-        self.length_scale = np.cbrt(self.particulator.environment.mesh.dv)
 
     def __call__(self):
         """for stochastic removal, see, e.g., the naive scheme in Algorithm 1 in
         [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)"""
+        self.length_scale = np.cbrt(self.particulator.environment.mesh.dv)
         self.particulator.sedimentation_removal(
             stochastic_sedimentation_removal=self.stochastic_sedimentation_removal,
             length_scale=self.length_scale,

From a6888159e718ff25f52868c1be3fcb73767e8d39 Mon Sep 17 00:00:00 2001
From: Sylwester Arabas <sylwester.arabas@agh.edu.pl>
Date: Mon, 1 Jun 2026 19:08:20 +0200
Subject: [PATCH 18/18] remove self.length_scale as it is not needed anymore

---
 PySDM/dynamics/sedimentation_removal_0d.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/PySDM/dynamics/sedimentation_removal_0d.py b/PySDM/dynamics/sedimentation_removal_0d.py
index 67dfbf8527..e26b18d8bf 100644
--- a/PySDM/dynamics/sedimentation_removal_0d.py
+++ b/PySDM/dynamics/sedimentation_removal_0d.py
@@ -10,7 +10,6 @@ class SedimentationRemoval0D:
     def __init__(self, *, stochastic_sedimentation_removal: Optional[bool] = True):
         """stochastic or deterministic removal"""
         self.stochastic_sedimentation_removal = stochastic_sedimentation_removal
-        self.length_scale = None
         self.particulator = None
 
     def register(self, builder):
@@ -21,8 +20,8 @@ def register(self, builder):
     def __call__(self):
         """for stochastic removal, see, e.g., the naive scheme in Algorithm 1 in
         [Curtis et al. 2016](https://doi.org/10.1016/j.jcp.2016.06.029)"""
-        self.length_scale = np.cbrt(self.particulator.environment.mesh.dv)
+
         self.particulator.sedimentation_removal(
             stochastic_sedimentation_removal=self.stochastic_sedimentation_removal,
-            length_scale=self.length_scale,
+            length_scale=np.cbrt(self.particulator.environment.mesh.dv),
         )