diff --git a/baybe/kernels/__init__.py b/baybe/kernels/__init__.py
index 9323a2b631..4e50c5728c 100644
--- a/baybe/kernels/__init__.py
+++ b/baybe/kernels/__init__.py
@@ -5,6 +5,7 @@
 """
 
 from baybe.kernels.basic import (
+    IndexKernel,
     LinearKernel,
     MaternKernel,
     PeriodicKernel,
@@ -18,6 +19,7 @@
 
 __all__ = [
     "AdditiveKernel",
+    "IndexKernel",
     "LinearKernel",
     "MaternKernel",
     "PeriodicKernel",
diff --git a/baybe/parameters/fidelity.py b/baybe/parameters/fidelity.py
index 15235c9c78..14ecf1a51a 100644
--- a/baybe/parameters/fidelity.py
+++ b/baybe/parameters/fidelity.py
@@ -21,8 +21,8 @@
     validate_is_finite,
     validate_unique_values,
 )
+from baybe.settings import active_settings
 from baybe.utils.conversion import nonstring_to_tuple
-from baybe.utils.numerical import DTypeFloatNumpy
 
 
 def _convert_zeta(
@@ -107,7 +107,9 @@ def values(self) -> tuple[str | bool, ...]:
     @cached_property
     def comp_df(self) -> pd.DataFrame:
         return pd.DataFrame(
-            range(len(self.values)), dtype=DTypeFloatNumpy, columns=[self.name]
+            range(len(self.values)),
+            dtype=active_settings.DTypeFloatNumpy,
+            columns=[self.name],
         )
 
 
@@ -159,5 +161,7 @@ def values(self) -> tuple[float, ...]:
     @cached_property
     def comp_df(self) -> pd.DataFrame:
         return pd.DataFrame(
-            {self.name: self.values}, index=self.values, dtype=DTypeFloatNumpy
+            {self.name: self.values},
+            index=self.values,
+            dtype=active_settings.DTypeFloatNumpy,
         )
diff --git a/baybe/recommenders/pure/bayesian/base.py b/baybe/recommenders/pure/bayesian/base.py
index 4ac5c1eed2..8da8607f6f 100644
--- a/baybe/recommenders/pure/bayesian/base.py
+++ b/baybe/recommenders/pure/bayesian/base.py
@@ -44,6 +44,13 @@ def _autoreplicate(surrogate: SurrogateProtocol, /) -> SurrogateProtocol:
 class BayesianRecommender(PureRecommender, ABC):
     """An abstract class for Bayesian Recommenders."""
 
+    # TODO: Factory defaults the surrogate to a GaussianProcessesSurrogate always.
+    #   Surrogate and kernel defaults should be different for searchspaces with
+    #   CategoricalFidelityParameter or NumericalDiscreteFidelityParameter.
+    #   This can be achieved without the user having to specify the surroagte model,
+    #   e.g., by
+    #   * using a dispatcher factory which decides surrogate model on fit time
+    #   * having a "_setup_surrogate" method similar to the acquisition function logic
     _surrogate_model: SurrogateProtocol = field(
         alias="surrogate_model",
         factory=GaussianProcessSurrogate,
diff --git a/baybe/searchspace/core.py b/baybe/searchspace/core.py
index 5510af704f..c2a47cbdd0 100644
--- a/baybe/searchspace/core.py
+++ b/baybe/searchspace/core.py
@@ -15,6 +15,10 @@
 from baybe.constraints.base import Constraint
 from baybe.parameters import TaskParameter
 from baybe.parameters.base import Parameter
+from baybe.parameters.fidelity import (
+    CategoricalFidelityParameter,
+    NumericalDiscreteFidelityParameter,
+)
 from baybe.searchspace.continuous import SubspaceContinuous
 from baybe.searchspace.discrete import (
     MemorySize,
@@ -48,6 +52,29 @@ class SearchSpaceType(Enum):
     """Flag for hybrid search spaces resp. compatibility with hybrid search spaces."""
 
 
+class SearchSpaceTaskType(Enum):
+    """Enum class for different types of task and/or fidelity subspaces."""
+
+    SINGLETASK = "SINGLETASK"
+    """Flag for search spaces with no task parameters."""
+
+    CATEGORICALTASK = "CATEGORICALTASK"
+    """Flag for search spaces with a categorical task parameter."""
+
+    NUMERICALFIDELITY = "NUMERICALFIDELITY"
+    """Flag for search spaces with a discrete numerical (ordered) fidelity parameter."""
+
+    CATEGORICALFIDELITY = "CATEGORICALFIDELITY"
+    """Flag for search spaces with a categorical (unordered) fidelity parameter."""
+
+    # TODO: Distinguish between multiple task parameter and mixed task parameter types.
+    # In future versions, multiple task/fidelity parameters may be allowed. For now,
+    # they are disallowed, whether the task-like parameters are different or the same
+    # class.
+    MULTIPLETASKPARAMETER = "MULTIPLETASKPARAMETER"
+    """Flag for search spaces with mixed task and fidelity parameters."""
+
+
 @define
 class SearchSpace(SerialMixin):
     """Class for managing the overall search space.
@@ -275,6 +302,24 @@ def task_idx(self) -> int | None:
         #   --> Fix this when refactoring the data
         return cast(int, self.discrete.comp_rep.columns.get_loc(task_param.name))
 
+    @property
+    def fidelity_idx(self) -> int | None:
+        """The column index of the task parameter in computational representation."""
+        try:
+            # See TODO [16932] and TODO [11611]
+            fidelity_param = next(
+                p
+                for p in self.parameters
+                if isinstance(
+                    p,
+                    (CategoricalFidelityParameter, NumericalDiscreteFidelityParameter),
+                )
+            )
+        except StopIteration:
+            return None
+
+        return cast(int, self.discrete.comp_rep.columns.get_loc(fidelity_param.name))
+
     @property
     def n_tasks(self) -> int:
         """The number of tasks encoded in the search space."""
@@ -287,6 +332,105 @@ def n_tasks(self) -> int:
             return 1
         return len(task_param.values)
 
+    @property
+    def n_fidelities(self) -> int:
+        """The number of tasks encoded in the search space."""
+        # See TODO [16932]
+        try:
+            fidelity_param = next(
+                p
+                for p in self.parameters
+                if isinstance(
+                    p,
+                    (CategoricalFidelityParameter, NumericalDiscreteFidelityParameter),
+                )
+            )
+            return len(fidelity_param.values)
+
+        # When there are no fidelity parameters, we effectively have a single fidelity
+        except StopIteration:
+            return 1
+
+    @property
+    def n_task_dimensions(self) -> int:
+        """The number of task dimensions."""
+        try:
+            # See TODO [16932]
+            fidelity_param = next(
+                p for p in self.parameters if isinstance(p, (TaskParameter,))
+            )
+        except StopIteration:
+            fidelity_param = None
+
+        return 1 if fidelity_param is not None else 0
+
+    @property
+    def n_fidelity_dimensions(self) -> int:
+        """The number of fidelity dimensions."""
+        try:
+            # See TODO [16932]
+            fidelity_param = next(
+                p
+                for p in self.parameters
+                if isinstance(
+                    p,
+                    (CategoricalFidelityParameter, NumericalDiscreteFidelityParameter),
+                )
+            )
+        except StopIteration:
+            fidelity_param = None
+
+        return 1 if fidelity_param is not None else 0
+
+    @property
+    def task_type(self) -> SearchSpaceTaskType:
+        """Return the task type of the search space.
+
+        Raises:
+            ValueError: If searchspace contains more than one task/fidelity parameter.
+            ValueError: An unrecognised fidelity parameter type is in SearchSpace.
+        """
+        task_like_parameters = (
+            TaskParameter,
+            CategoricalFidelityParameter,
+            NumericalDiscreteFidelityParameter,
+        )
+
+        n_task_like_parameters = sum(
+            isinstance(p, (task_like_parameters)) for p in self.parameters
+        )
+
+        if n_task_like_parameters == 0:
+            return SearchSpaceTaskType.SINGLETASK
+        elif n_task_like_parameters > 1:
+            # TODO: commute this validation further downstream.
+            # In case of user-defined custom models which allow for multiple task
+            # parameters, this should be later in recommender logic.
+            #  * Should this be an IncompatibilityError?
+            raise ValueError(
+                "SearchSpace must not contain more than one task/fidelity parameter."
+            )
+            return SearchSpaceTaskType.MULTIPLETASKPARAMETER
+
+        if self.n_task_dimensions == 1:
+            return SearchSpaceTaskType.CATEGORICALTASK
+
+        if self.n_fidelity_dimensions == 1:
+            n_categorical_fidelity_dims = sum(
+                isinstance(p, CategoricalFidelityParameter) for p in self.parameters
+            )
+            if n_categorical_fidelity_dims == 1:
+                return SearchSpaceTaskType.CATEGORICALFIDELITY
+
+            n_numerical_disc_fidelity_dims = sum(
+                isinstance(p, NumericalDiscreteFidelityParameter)
+                for p in self.parameters
+            )
+            if n_numerical_disc_fidelity_dims == 1:
+                return SearchSpaceTaskType.NUMERICALFIDELITY
+
+        raise RuntimeError("This line should be impossible to reach.")
+
     def get_comp_rep_parameter_indices(self, name: str, /) -> tuple[int, ...]:
         """Find a parameter's column indices in the computational representation.
 
diff --git a/baybe/surrogates/bandit.py b/baybe/surrogates/bandit.py
index ad6563cc43..3437e2b945 100644
--- a/baybe/surrogates/bandit.py
+++ b/baybe/surrogates/bandit.py
@@ -32,6 +32,9 @@ class BetaBernoulliMultiArmedBanditSurrogate(Surrogate):
     supports_transfer_learning: ClassVar[bool] = False
     # See base class.
 
+    supports_multi_fidelity: ClassVar[bool] = False
+    # See base class.
+
     prior: BetaPrior = field(factory=lambda: BetaPrior(1, 1))
     """The beta prior for the win rates of the bandit arms. Uniform by default."""
 
diff --git a/baybe/surrogates/base.py b/baybe/surrogates/base.py
index 205e32f703..244b0320e5 100644
--- a/baybe/surrogates/base.py
+++ b/baybe/surrogates/base.py
@@ -86,6 +86,10 @@ class Surrogate(ABC, SurrogateProtocol, SerialMixin):
     """Class variable encoding whether or not the surrogate supports transfer
     learning."""
 
+    supports_multi_fidelity: ClassVar[bool]
+    """Class variable encoding whether or not the surrogate supports multi fidelity
+    Bayesian optimization."""
+
     supports_multi_output: ClassVar[bool] = False
     """Class variable encoding whether or not the surrogate is multi-output
     compatible."""
@@ -428,6 +432,14 @@ def fit(
                 f"support transfer learning."
             )
 
+        # Check if multi fidelity capabilities are needed
+        if (searchspace.n_fidelities > 1) and (not self.supports_multi_fidelity):
+            raise ValueError(
+                f"The search space contains fidelity parameters but the selected "
+                f"surrogate model type ({self.__class__.__name__}) does not "
+                f"support multi fidelity Bayesian optimisation."
+            )
+
         # Block partial measurements
         handle_missing_values(measurements, [t.name for t in objective.targets])
 
@@ -472,6 +484,11 @@ def __str__(self) -> str:
                 self.supports_transfer_learning,
                 single_line=True,
             ),
+            to_string(
+                "Supports Multi Fidelity",
+                self.supports_multi_fidelity,
+                single_line=True,
+            ),
         ]
         return to_string(self.__class__.__name__, *fields)
 
diff --git a/baybe/surrogates/custom.py b/baybe/surrogates/custom.py
index 79c4c6ea1c..2b65b08a5f 100644
--- a/baybe/surrogates/custom.py
+++ b/baybe/surrogates/custom.py
@@ -70,6 +70,9 @@ class CustomONNXSurrogate(IndependentGaussianSurrogate):
     supports_transfer_learning: ClassVar[bool] = False
     # See base class.
 
+    supports_multi_fidelity: ClassVar[bool] = False
+    # See base class.
+
     onnx_input_name: str = field(validator=validators.instance_of(str))
     """The input name used for constructing the ONNX str."""
 
diff --git a/baybe/surrogates/gaussian_process/components/kernel.py b/baybe/surrogates/gaussian_process/components/kernel.py
index 1f038345ff..eb70ded893 100644
--- a/baybe/surrogates/gaussian_process/components/kernel.py
+++ b/baybe/surrogates/gaussian_process/components/kernel.py
@@ -10,6 +10,7 @@
 from baybe.kernels.base import Kernel
 from baybe.kernels.composite import ProductKernel
 from baybe.parameters.categorical import TaskParameter
+from baybe.parameters.fidelity import CategoricalFidelityParameter
 from baybe.parameters.selector import (
     ParameterSelectorProtocol,
     TypeSelector,
@@ -79,7 +80,15 @@ def _default_base_kernel_factory(self) -> KernelFactoryProtocol:
             BayBENumericalKernelFactory,
         )
 
-        return BayBENumericalKernelFactory(TypeSelector((TaskParameter,), exclude=True))
+        return BayBENumericalKernelFactory(
+            TypeSelector(
+                (
+                    TaskParameter,
+                    CategoricalFidelityParameter,
+                ),
+                exclude=True,
+            )
+        )
 
     @task_kernel_factory.default
     def _default_task_kernel_factory(self) -> KernelFactoryProtocol:
@@ -87,7 +96,14 @@ def _default_task_kernel_factory(self) -> KernelFactoryProtocol:
             BayBETaskKernelFactory,
         )
 
-        return BayBETaskKernelFactory(TypeSelector((TaskParameter,)))
+        return BayBETaskKernelFactory(
+            TypeSelector(
+                (
+                    TaskParameter,
+                    CategoricalFidelityParameter,
+                )
+            )
+        )
 
     @override
     def __call__(
diff --git a/baybe/surrogates/gaussian_process/core.py b/baybe/surrogates/gaussian_process/core.py
index 617a4a247c..2c70a09151 100644
--- a/baybe/surrogates/gaussian_process/core.py
+++ b/baybe/surrogates/gaussian_process/core.py
@@ -52,6 +52,8 @@
     from torch import Tensor
 
 
+# TODO Jordan MHS: _ModelContext is used by fidelity surrogate models now so may deserve
+# its own file.
 @define
 class _ModelContext:
     """Model context for :class:`GaussianProcessSurrogate`."""
@@ -80,6 +82,27 @@ def n_tasks(self) -> int:
         """The number of tasks."""
         return self.searchspace.n_tasks
 
+    @property
+    def n_fidelity_dimensions(self) -> int:
+        """The number of fidelity dimensions."""
+        # Possible TODO: Generalize to multiple fidelity dimensions
+        return 1 if self.searchspace.fidelity_idx is not None else 0
+
+    @property
+    def is_multi_fidelity(self) -> bool:
+        """Are there any fidelity dimensions?"""
+        return self.n_fidelity_dimensions > 0
+
+    @property
+    def fidelity_idx(self) -> int | None:
+        """The computational column index of the task parameter, if available."""
+        return self.searchspace.fidelity_idx
+
+    @property
+    def n_fidelities(self) -> int:
+        """The number of fidelities."""
+        return self.searchspace.n_fidelities
+
     @property
     def parameter_bounds(self) -> Tensor:
         """Get the search space parameter bounds in BoTorch Format."""
@@ -93,7 +116,7 @@ def numerical_indices(self) -> list[int]:
         return [
             i
             for i in range(len(self.searchspace.comp_rep_columns))
-            if i != self.task_idx
+            if i not in (self.task_idx, self.fidelity_idx)
         ]
 
 
diff --git a/baybe/surrogates/gaussian_process/multi_fidelity.py b/baybe/surrogates/gaussian_process/multi_fidelity.py
new file mode 100644
index 0000000000..2cdca81801
--- /dev/null
+++ b/baybe/surrogates/gaussian_process/multi_fidelity.py
@@ -0,0 +1,114 @@
+"""Multi-fidelity Gaussian process surrogates."""
+
+from __future__ import annotations
+
+import gc
+from typing import TYPE_CHECKING, ClassVar
+
+from attrs import define, field
+from typing_extensions import override
+
+from baybe.parameters.base import Parameter
+from baybe.surrogates.base import Surrogate
+from baybe.surrogates.gaussian_process.core import (
+    _ModelContext,
+)
+from baybe.surrogates.gaussian_process.presets.core import (
+    GaussianProcessPreset,
+    make_gp_from_preset,
+)
+
+if TYPE_CHECKING:
+    from botorch.models.gpytorch import GPyTorchModel
+    from botorch.models.transforms.input import InputTransform
+    from botorch.models.transforms.outcome import OutcomeTransform
+    from botorch.posteriors import Posterior
+    from torch import Tensor
+
+
+@define
+class GaussianProcessSurrogateSTMF(Surrogate):
+    """Botorch's single task multi fidelity Gaussian process."""
+
+    supports_transfer_learning: ClassVar[bool] = False
+    # See base class.
+
+    supports_multi_fidelity: ClassVar[bool] = True
+    # See base class.
+
+    # TODO: type should be Optional[botorch.models.SingleTaskGP] but is currently
+    #   omitted due to: https://github.com/python-attrs/cattrs/issues/531
+    _model = field(init=False, default=None, eq=False)
+    """The actual model."""
+
+    @staticmethod
+    def from_preset(preset: GaussianProcessPreset) -> Surrogate:
+        """Create a Gaussian process surrogate from one of the defined presets."""
+        return make_gp_from_preset(preset)
+
+    @override
+    def to_botorch(self) -> GPyTorchModel:
+        return self._model
+
+    @override
+    @staticmethod
+    def _make_parameter_scaler_factory(
+        parameter: Parameter,
+    ) -> type[InputTransform] | None:
+        # For GPs, we let botorch handle the scaling. See [Scaling Workaround] above.
+        return None
+
+    @override
+    @staticmethod
+    def _make_target_scaler_factory() -> type[OutcomeTransform] | None:
+        # For GPs, we let botorch handle the scaling. See [Scaling Workaround] above.
+        return None
+
+    @override
+    def _posterior(self, candidates_comp_scaled: Tensor, /) -> Posterior:
+        return self._model.posterior(candidates_comp_scaled)
+
+    @override
+    def _fit(self, train_x: Tensor, train_y: Tensor) -> None:
+        import botorch
+        import gpytorch
+
+        assert self._searchspace is not None
+
+        context = _ModelContext(self._searchspace)
+
+        assert context.is_multi_fidelity, (
+            "GaussianProcessSurrogateSTMF can only "
+            "be fit on multi fidelity searchspaces."
+        )
+
+        # For GPs, we let botorch handle the scaling. See [Scaling Workaround] above.
+        input_transform = botorch.models.transforms.Normalize(  # type: ignore[attr-defined]
+            train_x.shape[-1],
+            bounds=context.parameter_bounds,
+            indices=context.numerical_indices,
+        )
+        outcome_transform = botorch.models.transforms.Standardize(train_y.shape[-1])  # type: ignore[attr-defined]
+
+        # construct and fit the Gaussian process
+        self._model = botorch.models.SingleTaskMultiFidelityGP(
+            train_x,
+            train_y,
+            input_transform=input_transform,
+            outcome_transform=outcome_transform,
+            data_fidelities=None
+            if context.fidelity_idx is None
+            else (context.fidelity_idx,),
+        )
+
+        mll = gpytorch.ExactMarginalLogLikelihood(self._model.likelihood, self._model)
+
+        botorch.fit.fit_gpytorch_mll(mll)
+
+    @override
+    def __str__(self) -> str:
+        return "SingleTaskMultiFidelityGP with Botorch defaults."
+
+
+# Collect leftover original slotted classes processed by `attrs.define`
+gc.collect()
diff --git a/baybe/surrogates/gaussian_process/presets/core.py b/baybe/surrogates/gaussian_process/presets/core.py
index ad77df0b4d..5f659c0e01 100644
--- a/baybe/surrogates/gaussian_process/presets/core.py
+++ b/baybe/surrogates/gaussian_process/presets/core.py
@@ -3,6 +3,10 @@
 from __future__ import annotations
 
 from enum import Enum
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from baybe.surrogates.base import Surrogate
 
 
 class GaussianProcessPreset(Enum):
@@ -16,3 +20,20 @@ class GaussianProcessPreset(Enum):
 
     EDBO_SMOOTHED = "EDBO_SMOOTHED"
     """A smoothed version of the EDBO settings."""
+
+    BOTORCH_STMF = "BOTORCH_STMF"
+    """Recreates the default settings of the BOTORCH SingleTaskMultiFidelityGP."""
+
+
+def make_gp_from_preset(preset: GaussianProcessPreset) -> Surrogate:
+    """Create a :class:`GaussianProcessSurrogate` from a :class:`GaussianProcessPreset."""  # noqa: E501
+    from baybe.surrogates.gaussian_process.multi_fidelity import (
+        GaussianProcessSurrogateSTMF,
+    )
+
+    if preset is GaussianProcessPreset.BOTORCH_STMF:
+        return GaussianProcessSurrogateSTMF()
+
+    raise ValueError(
+        f"Unknown '{GaussianProcessPreset.__name__}' with name '{preset.name}'."
+    )
diff --git a/baybe/surrogates/linear.py b/baybe/surrogates/linear.py
index ba1fed5b2c..94746a16b9 100644
--- a/baybe/surrogates/linear.py
+++ b/baybe/surrogates/linear.py
@@ -44,6 +44,9 @@ class BayesianLinearSurrogate(IndependentGaussianSurrogate):
     supports_transfer_learning: ClassVar[bool] = False
     # See base class.
 
+    supports_multi_fidelity: ClassVar[bool] = False
+    # See base class.
+
     model_params: _ARDRegressionParams = field(
         factory=dict,
         converter=dict,
diff --git a/baybe/surrogates/naive.py b/baybe/surrogates/naive.py
index 3912c6b128..b407b48f08 100644
--- a/baybe/surrogates/naive.py
+++ b/baybe/surrogates/naive.py
@@ -26,6 +26,9 @@ class MeanPredictionSurrogate(IndependentGaussianSurrogate):
     supports_transfer_learning: ClassVar[bool] = False
     # See base class.
 
+    supports_multi_fidelity: ClassVar[bool] = False
+    # See base class.
+
     _model: float | None = field(init=False, default=None, eq=False)
     """The estimated posterior mean value of the training targets."""
 
diff --git a/baybe/surrogates/ngboost.py b/baybe/surrogates/ngboost.py
index 6cb1cee135..f05a9ebc0d 100644
--- a/baybe/surrogates/ngboost.py
+++ b/baybe/surrogates/ngboost.py
@@ -49,6 +49,9 @@ class NGBoostSurrogate(IndependentGaussianSurrogate):
     supports_transfer_learning: ClassVar[bool] = False
     # See base class.
 
+    supports_multi_fidelity: ClassVar[bool] = False
+    # See base class.
+
     _default_model_params: ClassVar[dict] = {"n_estimators": 25, "verbose": False}
     """Class variable encoding the default model parameters."""
 
diff --git a/baybe/surrogates/random_forest.py b/baybe/surrogates/random_forest.py
index 91ad599bb1..6ee1f7ed70 100644
--- a/baybe/surrogates/random_forest.py
+++ b/baybe/surrogates/random_forest.py
@@ -64,6 +64,9 @@ class RandomForestSurrogate(Surrogate):
     supports_transfer_learning: ClassVar[bool] = False
     # See base class.
 
+    supports_multi_fidelity: ClassVar[bool] = False
+    # See base class.
+
     model_params: _RandomForestRegressorParams = field(
         factory=dict,
         converter=dict,