diff --git a/CITATION.cff b/CITATION.cff
index 3dcf0f2..87b809b 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -8,7 +8,7 @@ authors:
 - family-names: "Weissgraeber"
   given-names: "Philipp"
   orcid: "https://orcid.org/0000-0001-8320-8672"
-version: 3.1.2
+version: 3.1.3
 date-released: 2021-12-30
 identifiers:
   - description: Collection of archived snapshots of all versions of WEAC
diff --git a/demo/demo.ipynb b/demo/demo.ipynb
index 38007da..df82164 100644
--- a/demo/demo.ipynb
+++ b/demo/demo.ipynb
@@ -13,7 +13,7 @@
    "id": "695bafcb",
    "metadata": {},
    "source": [
-    "Note that instructions in this notebook refer to **release v3.1.2.** Please make sure you are running the latest version of weac using\n",
+    "Note that instructions in this notebook refer to **release v3.1.3.** Please make sure you are running the latest version of weac using\n",
     "\n",
     "```bash\n",
     "pip install -U weac\n",
diff --git a/pyproject.toml b/pyproject.toml
index cf35b60..3e242e7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "weac"
-version = "3.1.2"
+version = "3.1.3"
 authors = [{ name = "2phi GbR", email = "mail@2phi.de" }]
 description = "Weak layer anticrack nucleation model"
 readme = "README.md"
@@ -123,7 +123,7 @@ ignore = [
 ]
 
 [tool.bumpversion]
-current_version = "3.1.2"
+current_version = "3.1.3"
 
 [[tool.bumpversion.files]]
 filename = "pyproject.toml"
diff --git a/src/weac/__init__.py b/src/weac/__init__.py
index 58ccd86..a9798b2 100644
--- a/src/weac/__init__.py
+++ b/src/weac/__init__.py
@@ -2,4 +2,4 @@
 WEAC - Weak Layer Anticrack Nucleation Model
 """
 
-__version__ = "3.1.2"
+__version__ = "3.1.3"
diff --git a/src/weac/analysis/criteria_evaluator.py b/src/weac/analysis/criteria_evaluator.py
index 90de157..006bce9 100644
--- a/src/weac/analysis/criteria_evaluator.py
+++ b/src/weac/analysis/criteria_evaluator.py
@@ -24,7 +24,7 @@
     WeakLayer,
 )
 from weac.constants import RHO_ICE
-from weac.core.slab_touchdown import TouchdownMode
+from weac.components.scenario_config import TouchdownMode
 from weac.core.system_model import SystemModel
 
 logger = logging.getLogger(__name__)
@@ -741,6 +741,9 @@ def evaluate_SteadyState(
             cut_length=cut_distance,
         )
         system_copy.update_scenario(segments=segments, scenario_config=scenario_config)
+        # Force the requested mode to avoid floating-point precision issues
+        # when l_AB/l_BC are recalculated with different scenario parameters
+        system_copy.set_forced_touchdown_mode(mode)
         touchdown_distance = system_copy.slab_touchdown.touchdown_distance
         analyzer = Analyzer(system_copy, printing_enabled=print_call_stats)
         energy_release_rate, _, _ = analyzer.differential_ERR(unit="J/m^2")
diff --git a/src/weac/components/__init__.py b/src/weac/components/__init__.py
index 0c199bb..724514b 100644
--- a/src/weac/components/__init__.py
+++ b/src/weac/components/__init__.py
@@ -7,7 +7,7 @@
 from .layer import Layer, WeakLayer
 from .model_input import ModelInput
 from .segment import Segment
-from .scenario_config import ScenarioConfig, SystemType
+from .scenario_config import ScenarioConfig, SystemType, TouchdownMode
 
 __all__ = [
     "Config",
@@ -18,4 +18,5 @@
     "ScenarioConfig",
     "ModelInput",
     "SystemType",
+    "TouchdownMode",
 ]
diff --git a/src/weac/components/config.py b/src/weac/components/config.py
index 590c974..7ce849a 100644
--- a/src/weac/components/config.py
+++ b/src/weac/components/config.py
@@ -12,6 +12,8 @@
 
 from pydantic import BaseModel, Field
 
+from weac.components.scenario_config import TouchdownMode
+
 
 class Config(BaseModel):
     """
@@ -21,11 +23,19 @@ class Config(BaseModel):
     ----------
     touchdown : bool
         Whether slab touchdown on the collapsed weak layer is considered.
+    forced_touchdown_mode : TouchdownMode | None
+        If set, forces the touchdown mode instead of calculating it from l_AB/l_BC.
+        This avoids floating-point precision issues when the mode boundary values
+        are recalculated with different scenario parameters.
     """
 
     touchdown: bool = Field(
         default=False, description="Whether to include slab touchdown in the analysis"
     )
+    forced_touchdown_mode: TouchdownMode | None = Field(
+        default=None,
+        description="Force a specific touchdown mode instead of auto-calculating",
+    )
 
 
 if __name__ == "__main__":
diff --git a/src/weac/components/scenario_config.py b/src/weac/components/scenario_config.py
index 17fccaa..0b6b715 100644
--- a/src/weac/components/scenario_config.py
+++ b/src/weac/components/scenario_config.py
@@ -11,6 +11,8 @@
     "skier", "skiers", "pst-", "-pst", "rot", "trans", "vpst-", "-vpst"
 ]
 
+TouchdownMode = Literal["A_free_hanging", "B_point_contact", "C_in_contact"]
+
 
 class ScenarioConfig(BaseModel):
     """
diff --git a/src/weac/core/slab_touchdown.py b/src/weac/core/slab_touchdown.py
index 36fe09f..66abb4f 100644
--- a/src/weac/core/slab_touchdown.py
+++ b/src/weac/core/slab_touchdown.py
@@ -9,7 +9,7 @@
 from scipy.optimize import brentq
 
 from weac.components.layer import WeakLayer
-from weac.components.scenario_config import ScenarioConfig
+from weac.components.scenario_config import ScenarioConfig, TouchdownMode
 from weac.components.segment import Segment
 from weac.constants import STIFFNESS_COLLAPSE_FACTOR
 from weac.core.eigensystem import Eigensystem
@@ -20,9 +20,6 @@
 logger = logging.getLogger(__name__)
 
 
-TouchdownMode = Literal["A_free_hanging", "B_point_contact", "C_in_contact"]
-
-
 class SlabTouchdown:  # pylint: disable=too-many-instance-attributes,too-few-public-methods
     """
     Handling the touchdown situation in a PST.
@@ -52,6 +49,9 @@ class SlabTouchdown:  # pylint: disable=too-many-instance-attributes,too-few-pub
     -----------
     scenario: `Scenario`
     eigensystem: `Eigensystem`
+    forced_mode: `TouchdownMode | None`
+        If provided, forces this touchdown mode instead of calculating from l_AB/l_BC.
+        This avoids floating-point precision issues when scenario parameters change.
 
     Attributes:
     -----------
@@ -81,9 +81,15 @@ class SlabTouchdown:  # pylint: disable=too-many-instance-attributes,too-few-pub
     touchdown_distance: float
     collapsed_weak_layer_kR: float | None = None
 
-    def __init__(self, scenario: Scenario, eigensystem: Eigensystem):
+    def __init__(
+        self,
+        scenario: Scenario,
+        eigensystem: Eigensystem,
+        forced_mode: TouchdownMode | None = None,
+    ):
         self.scenario = scenario
         self.eigensystem = eigensystem
+        self._forced_mode = forced_mode
 
         # Create a new scenario config with phi=0 (flat slab) while preserving other settings
         self.flat_config = ScenarioConfig(
@@ -114,14 +120,25 @@ def _calc_touchdown_mode(self):
             self.l_BC = self._calc_l_BC()
         except ValueError:
             self.l_BC = self.scenario.L
-        # Assign stage
-        touchdown_mode = "A_free_hanging"
-        if self.scenario.cut_length <= self.l_AB:
+
+        # Assign stage - use forced mode if provided, otherwise calculate from thresholds
+        if self._forced_mode is not None:
+            touchdown_mode = self._forced_mode
+            logger.debug(
+                "Using forced touchdown mode: %s (l_AB=%.2f, l_BC=%.2f, cut_length=%.2f)",
+                touchdown_mode,
+                self.l_AB,
+                self.l_BC,
+                self.scenario.cut_length,
+            )
+        else:
             touchdown_mode = "A_free_hanging"
-        elif self.l_AB < self.scenario.cut_length <= self.l_BC:
-            touchdown_mode = "B_point_contact"
-        elif self.l_BC < self.scenario.cut_length:
-            touchdown_mode = "C_in_contact"
+            if self.scenario.cut_length <= self.l_AB:
+                touchdown_mode = "A_free_hanging"
+            elif self.l_AB < self.scenario.cut_length <= self.l_BC:
+                touchdown_mode = "B_point_contact"
+            elif self.l_BC < self.scenario.cut_length:
+                touchdown_mode = "C_in_contact"
         self.touchdown_mode = touchdown_mode
 
     def _calc_touchdown_distance(self):
diff --git a/src/weac/core/system_model.py b/src/weac/core/system_model.py
index c9a8426..4ce5712 100644
--- a/src/weac/core/system_model.py
+++ b/src/weac/core/system_model.py
@@ -28,6 +28,7 @@
 from weac.core.field_quantities import FieldQuantities
 from weac.core.scenario import Scenario
 from weac.core.slab import Slab
+from weac.components.scenario_config import TouchdownMode
 from weac.core.slab_touchdown import SlabTouchdown
 from weac.core.unknown_constants_solver import UnknownConstantsSolver
 
@@ -164,7 +165,9 @@ def slab_touchdown(self) -> SlabTouchdown | None:
         if self.config.touchdown:
             logger.info("Solving for Slab Touchdown")
             slab_touchdown = SlabTouchdown(
-                scenario=self.scenario, eigensystem=self.eigensystem
+                scenario=self.scenario,
+                eigensystem=self.eigensystem,
+                forced_mode=self.config.forced_touchdown_mode,
             )
             logger.info(
                 "Original cut_length: %s, touchdown_distance: %s",
@@ -345,6 +348,26 @@ def toggle_touchdown(self, touchdown: bool):
             self._invalidate_slab_touchdown()
             self._invalidate_constants()
 
+    def set_forced_touchdown_mode(self, mode: TouchdownMode | None):
+        """
+        Set the forced touchdown mode.
+
+        When set, the touchdown mode calculation will use this mode instead of
+        calculating it from l_AB/l_BC thresholds. This avoids floating-point
+        precision issues when the mode boundaries are recalculated with different
+        scenario parameters.
+
+        Parameters
+        ----------
+        mode : TouchdownMode | None
+            The mode to force ("A_free_hanging", "B_point_contact", "C_in_contact"),
+            or None to use automatic calculation.
+        """
+        if self.config.forced_touchdown_mode != mode:
+            self.config.forced_touchdown_mode = mode
+            self._invalidate_slab_touchdown()
+            self._invalidate_constants()
+
     def _invalidate_eigensystem(self):
         """Invalidate the eigensystem."""
         self.__dict__.pop("eigensystem", None)
diff --git a/uv.lock b/uv.lock
index 2f0ab71..666984d 100644
--- a/uv.lock
+++ b/uv.lock
@@ -2809,7 +2809,7 @@ wheels = [
 
 [[package]]
 name = "weac"
-version = "3.1.1"
+version = "3.1.3"
 source = { editable = "." }
 dependencies = [
     { name = "matplotlib" },