diff --git a/pyproject.toml b/pyproject.toml
index 69db440a..fba38d13 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -55,7 +55,7 @@ exclude = '''
 
 [tool.ruff]
 line-length = 88
-select = ["E", "F", "I", "Q", "W"]
+select = ["E", "F", "I", "Q", "W", "PT"]
 # line too long; Black will handle these
 ignore = ["E501"]
 
diff --git a/src/calliope/backend/backend_model.py b/src/calliope/backend/backend_model.py
index d9ff3532..97e1a4a7 100644
--- a/src/calliope/backend/backend_model.py
+++ b/src/calliope/backend/backend_model.py
@@ -7,7 +7,6 @@
 
 from __future__ import annotations
 
-import importlib
 import logging
 import time
 import typing
@@ -36,11 +35,9 @@
 from calliope.exceptions import warn as model_warn
 from calliope.io import load_config
 from calliope.util.schema import (
-    MATH_SCHEMA,
     MODEL_SCHEMA,
     extract_from_schema,
     update_then_validate_config,
-    validate_dict,
 )
 
 if TYPE_CHECKING:
@@ -211,7 +208,6 @@ def _check_inputs(self):
         )
 
     def _build(self) -> None:
-        self._add_run_mode_math()
         # The order of adding components matters!
         # 1. Variables, 2. Global Expressions, 3. Constraints, 4. Objectives
         for components in [
@@ -230,27 +226,6 @@ def _build(self) -> None:
                 )
             LOGGER.info(f"Optimisation Model | {components} | Generated.")
 
-    def _add_run_mode_math(self) -> None:
-        """If not given in the add_math list, override model math with run mode math"""
-
-        # FIXME: available modes should not be hardcoded here. They should come from a YAML schema.
-        mode = self.inputs.attrs["config"].build.mode
-        add_math = self.inputs.attrs["applied_additional_math"]
-        not_run_mode = {"plan", "operate", "spores"}.difference([mode])
-        run_mode_mismatch = not_run_mode.intersection(add_math)
-        if run_mode_mismatch:
-            exceptions.warn(
-                f"Running in {mode} mode, but run mode(s) {run_mode_mismatch} "
-                "math being loaded from file via the model configuration"
-            )
-
-        if mode != "plan" and mode not in add_math:
-            LOGGER.debug(f"Updating math formulation with {mode} mode math.")
-            filepath = importlib.resources.files("calliope") / "math" / f"{mode}.yaml"
-            self.inputs.math.union(AttrDict.from_yaml(filepath), allow_override=True)
-
-        validate_dict(self.inputs.math, MATH_SCHEMA, "math")
-
     def _add_component(
         self,
         name: str,
diff --git a/src/calliope/config/config_schema.yaml b/src/calliope/config/config_schema.yaml
index ca8749b0..e1ddc450 100644
--- a/src/calliope/config/config_schema.yaml
+++ b/src/calliope/config/config_schema.yaml
@@ -51,6 +51,13 @@ properties:
             type: string
             default: "ISO8601"
             description: Timestamp format of all time series data when read from file. "ISO8601" means "%Y-%m-%d %H:%M:%S".
+          base_math:
+            type: boolean
+            default: true
+            description: >-
+              Whether or not to use the base Calliope math and schema when building models.
+              If "true", the model will use the base Calliope math. Values in "add_math" will be merged into it.
+              If "false", users must fully define the model math via "add_math". Use with caution!
           add_math:
             type: array
             default: []
diff --git a/src/calliope/model.py b/src/calliope/model.py
index 68803b5a..a8e27256 100644
--- a/src/calliope/model.py
+++ b/src/calliope/model.py
@@ -188,9 +188,6 @@ def _init_from_model_def_dict(
         model_config.union(model_definition.pop("config"), allow_override=True)
 
         init_config = update_then_validate_config("init", model_config)
-        # We won't store `init` in `self.config`, so we pop it out now.
-        model_config.pop("init")
-
         if init_config["time_cluster"] is not None:
             init_config["time_cluster"] = relative_path(
                 self._model_def_path, init_config["time_cluster"]
@@ -223,9 +220,7 @@ def _init_from_model_def_dict(
             init_config, model_definition, data_sources, attributes, param_metadata
         )
         model_data_factory.build()
-
         self._model_data = model_data_factory.dataset
-        self._model_data.attrs["math"] = self._add_math(init_config["add_math"])
 
         log_time(
             LOGGER,
@@ -233,11 +228,17 @@ def _init_from_model_def_dict(
             "model_data_creation",
             comment="Model: preprocessing stage 2 (model_data)",
         )
+
+        math, applied_math = self._math_init_from_yaml(init_config)
+        self._model_data.attrs["math"] = math
+        self._model_data.attrs["applied_math"] = applied_math
+        validate_dict(self.math, MATH_SCHEMA, "math")
+
         log_time(
             LOGGER,
             self._timings,
             "model_preprocessing_complete",
-            comment="Model: preprocessing complete",
+            comment="Model: preprocessing complete (model_math)",
         )
 
     def _init_from_model_data(self, model_data: xr.Dataset) -> None:
@@ -273,27 +274,26 @@ def _init_from_model_data(self, model_data: xr.Dataset) -> None:
             comment="Model: loaded model_data",
         )
 
-    def _add_math(self, add_math: list) -> AttrDict:
-        """
-        Load the base math and optionally override with additional math from a list of references to math files.
+    def _math_init_from_yaml(self, init_config: dict) -> tuple[AttrDict, list]:
+        """Construct math by combining files specified in the initialisation configuration.
 
         Args:
-            add_math (list):
-                List of references to files containing mathematical formulations that will be merged with the base formulation.
+            init_config (dict): initialization configuration to apply to model math.
 
         Raises:
-            exceptions.ModelError:
-                Referenced pre-defined math files or user-defined math files must exist.
+            exceptions.ModelError: referenced pre-defined math files or user-defined math files must exist.
 
         Returns:
-            AttrDict: Dictionary of math (constraints, variables, objectives, and global expressions).
+            tuple[AttrDict, list]: initialized mathematical formulation and list of applied math files.
         """
-        math_dir = Path(calliope.__file__).parent / "math"
-        base_math = AttrDict.from_yaml(math_dir / "base.yaml")
+        applied_math = init_config["add_math"].copy()
+        if init_config["base_math"]:
+            applied_math.insert(0, "base")
 
+        math_dir = Path(calliope.__file__).parent / "math"
+        math = AttrDict()
         file_errors = []
-
-        for filename in add_math:
+        for filename in applied_math:
             if not f"{filename}".endswith((".yaml", ".yml")):
                 yaml_filepath = math_dir / f"{filename}.yaml"
             else:
@@ -302,26 +302,62 @@ def _add_math(self, add_math: list) -> AttrDict:
             if not yaml_filepath.is_file():
                 file_errors.append(filename)
                 continue
-            else:
-                override_dict = AttrDict.from_yaml(yaml_filepath)
+            math.union(AttrDict.from_yaml(yaml_filepath), allow_override=True)
 
-            base_math.union(override_dict, allow_override=True)
         if file_errors:
             raise exceptions.ModelError(
                 f"Attempted to load additional math that does not exist: {file_errors}"
             )
-        self._model_data.attrs["applied_additional_math"] = add_math
-        return base_math
+        return (math, applied_math)
+
+    def _math_update_with_mode(
+        self, build_config: AttrDict
+    ) -> tuple[Union[AttrDict, None], list]:
+        """Evaluate if the model's math needs to be updated to enable the run mode.
+
+        The updated math and new list of applied math will be returned after some checks.
+        Will return an empty dictionary if no update is needed.
+
+        Args:
+            build_config (AttrDict): build configuration to apply to model math.
+
+        Returns:
+            tuple[Union[AttrDict, None], list]: math update and new list of applied files.
+        """
+        # FIXME: available modes should not be hardcoded here. They should come from a YAML schema.
+        mode = build_config["mode"]
+        applied_math = self._model_data.attrs["applied_math"].copy()
+        update_math = None
+
+        not_run_mode = {"plan", "operate", "spores"}.difference([mode])
+        run_mode_mismatch = not_run_mode.intersection(applied_math)
+        if run_mode_mismatch:
+            exceptions.warn(
+                f"Running in {mode} mode, but run mode(s) {run_mode_mismatch} "
+                "math being loaded from file via the model configuration"
+            )
+
+        if mode != "plan" and mode not in applied_math:
+            LOGGER.debug(f"Updating math formulation with {mode} mode math.")
+            filepath = Path(calliope.__file__).parent / "math" / f"{mode}.yaml"
+            mode_math = AttrDict.from_yaml(filepath)
+
+            update_math = self._model_data.attrs["math"].copy()
+            update_math.union(mode_math, allow_override=True)
+            applied_math.append(mode)
+
+        return (update_math, applied_math)
 
     def build(self, force: bool = False, **kwargs) -> None:
         """Build description of the optimisation problem in the chosen backend interface.
 
         Args:
-            force (bool, optional):
-                If ``force`` is True, any existing results will be overwritten.
-                Defaults to False.
-        """
+            force (bool, optional): If True, any existing results will be overwritten. Defaults to False.
 
+        Raises:
+            exceptions.ModelError: trying to re-build a model without enabling the force flag.
+            exceptions.ModelError: trying to run operate mode with incompatible added math.
+        """
         if self._is_built and not force:
             raise exceptions.ModelError(
                 "This model object already has a built optimisation problem. Use model.build(force=True) "
@@ -347,6 +383,14 @@ def build(self, force: bool = False, **kwargs) -> None:
             )
         else:
             input = self._model_data
+
+        # Build updates are valid, update configuration and math if necessary
+        self._model_data.attrs["config"]["build"] = updated_build_config
+        math_update, applied_math = self._math_update_with_mode(updated_build_config)
+        if math_update is not None:
+            self._model_data.attrs["math"] = math_update
+            self._model_data.attrs["applied_math"] = applied_math
+
         backend_name = updated_build_config["backend"]
         backend = self._BACKENDS[backend_name](input, **updated_build_config)
         backend._build()
diff --git a/tests/test_backend_pyomo.py b/tests/test_backend_pyomo.py
index 64683cd8..cb3bf004 100755
--- a/tests/test_backend_pyomo.py
+++ b/tests/test_backend_pyomo.py
@@ -1,6 +1,4 @@
-import importlib
 import logging
-from copy import deepcopy
 from itertools import product
 
 import calliope.exceptions as exceptions
@@ -9,8 +7,6 @@
 import pyomo.kernel as pmo
 import pytest  # noqa: F401
 import xarray as xr
-from calliope.attrdict import AttrDict
-from calliope.backend.pyomo_backend_model import PyomoBackendModel
 
 from .common.util import build_test_model as build_model
 from .common.util import check_error_or_warning, check_variable_exists
@@ -18,7 +14,7 @@
 
 @pytest.mark.xfail(reason="Not expecting operate mode to work at the moment")
 class TestChecks:
-    @pytest.mark.parametrize("on", (True, False))
+    @pytest.mark.parametrize("on", [True, False])
     def test_operate_cyclic_storage(self, on):
         """Cannot have cyclic storage in operate mode"""
         if on is True:
@@ -60,7 +56,7 @@ def test_loading_timeseries_operate_efficiencies(self, param):
         with pytest.warns(exceptions.ModelWarning):
             m.build()  # will fail to complete run if there's a problem
 
-    @pytest.mark.parametrize("constr", ("max", "equals"))
+    @pytest.mark.parametrize("constr", ["max", "equals"])
     def test_operate_flow_out_min_relative(self, constr):
         """If we depend on a finite flow_cap, we have to error on a user failing to define it"""
         m = build_model(
@@ -84,7 +80,7 @@ def test_operate_flow_out_min_relative(self, constr):
             error, ["Operate mode: User must define a finite flow_cap"]
         )
 
-    @pytest.mark.parametrize("constr", ("max", "equals"))
+    @pytest.mark.parametrize("constr", ["max", "equals"])
     def test_operate_flow_cap_source_unit(self, constr):
         """If we depend on a finite flow_cap, we have to error on a user failing to define it"""
         m = build_model(
@@ -112,7 +108,7 @@ def test_operate_flow_cap_source_unit(self, constr):
                 m.build()
 
     @pytest.mark.parametrize(
-        ["source_unit", "constr"],
+        ("source_unit", "constr"),
         list(product(("absolute", "per_cap", "per_area"), ("max", "equals"))),
     )
     def test_operate_source_unit_with_area_use(self, source_unit, constr):
@@ -237,7 +233,7 @@ def test_operate_storage(self, param):
             ],
         )
 
-    @pytest.mark.parametrize("on", (True, False))
+    @pytest.mark.parametrize("on", [True, False])
     def test_operate_source_cap_max(self, on):
         """Some constraints, if not defined, will throw a warning and possibly change values in model_data"""
 
@@ -264,7 +260,7 @@ def test_operate_source_cap_max(self, on):
             )
             assert m._model_data.source_cap.loc["a", "test_supply_plus"].item() == 1e6
 
-    @pytest.mark.parametrize("on", (True, False))
+    @pytest.mark.parametrize("on", [True, False])
     def test_operate_storage_initial(self, on):
         """Some constraints, if not defined, will throw a warning and possibly change values in model_data"""
 
@@ -512,8 +508,8 @@ def test_loc_techs_not_cost_var_constraint(self, simple_conversion):
         assert "cost_var" not in simple_conversion.backend.expressions
 
     @pytest.mark.parametrize(
-        "tech,scenario,cost",
-        (
+        ("tech", "scenario", "cost"),
+        [
             ("test_supply_elec", "simple_supply", "flow_out"),
             ("test_supply_elec", "simple_supply", "flow_in"),
             ("test_supply_plus", "simple_supply_and_supply_plus", "flow_in"),
@@ -521,7 +517,7 @@ def test_loc_techs_not_cost_var_constraint(self, simple_conversion):
             ("test_transmission_elec", "simple_supply", "flow_out"),
             ("test_conversion", "simple_conversion", "flow_in"),
             ("test_conversion_plus", "simple_conversion_plus", "flow_out"),
-        ),
+        ],
     )
     def test_loc_techs_cost_var_constraint(self, tech, scenario, cost):
         """
@@ -672,7 +668,7 @@ def test_loc_techs_storage_capacity_milp_constraint(self):
         assert "storage_capacity" not in m.backend.constraints
 
     @pytest.mark.parametrize(
-        "scenario,tech,override",
+        ("scenario", "tech", "override"),
         [
             i + (j,)
             for i in [
@@ -697,7 +693,7 @@ def test_loc_techs_flow_capacity_storage_constraint(self, scenario, tech, overri
         )
 
     @pytest.mark.filterwarnings("ignore:(?s).*Integer:calliope.exceptions.ModelWarning")
-    @pytest.mark.parametrize("override", (("max", "min")))
+    @pytest.mark.parametrize("override", (["max", "min"]))
     def test_loc_techs_flow_capacity_milp_storage_constraint(self, override):
         """
         i for i in sets.loc_techs_store if constraint_exists(model_run, i, 'constraints.flow_cap_per_storage_cap_max')
@@ -732,7 +728,7 @@ def test_no_loc_techs_flow_capacity_storage_constraint(self, caplog):
             ]
         )
 
-    @pytest.mark.parametrize("override", ((None, "max", "min")))
+    @pytest.mark.parametrize("override", ([None, "max", "min"]))
     def test_loc_techs_resource_capacity_constraint(self, override):
         """
         i for i in sets.loc_techs_finite_resource_supply_plus
@@ -939,8 +935,8 @@ def test_loc_techs_flow_capacity_constraint_warning_on_infinite_equals(self):
         override = {
             "nodes.a.techs.test_supply_elec.constraints.flow_cap_equals": np.inf
         }
+        m = build_model(override, "simple_supply,two_hours,investment_costs")
         with pytest.raises(exceptions.ModelError) as error:
-            m = build_model(override, "simple_supply,two_hours,investment_costs")
             m.build()
 
         assert check_error_or_warning(
@@ -995,7 +991,7 @@ def check_bounds(constraint):
             ).sel(techs="test_supply_elec")
         )
 
-    @pytest.mark.parametrize("bound", (("equals", "max")))
+    @pytest.mark.parametrize("bound", (["equals", "max"]))
     def test_techs_flow_capacity_systemwide_no_constraint(self, simple_supply, bound):
         assert "flow_capacity_systemwide" not in simple_supply.backend.constraints
 
@@ -1433,7 +1429,7 @@ def test_loc_techs_storage_capacity_min_purchase_milp_constraint(
 
     @pytest.mark.parametrize(
         ("scenario", "exists", "override_dict"),
-        (
+        [
             ("simple_supply", ("not", "not"), {}),
             ("supply_milp", ("not", "not"), {}),
             (
@@ -1442,7 +1438,7 @@ def test_loc_techs_storage_capacity_min_purchase_milp_constraint(
                 {"techs.test_supply_elec.costs.monetary.purchase": 1},
             ),
             ("supply_purchase", ("is", "not"), {}),
-        ),
+        ],
     )
     def test_loc_techs_update_costs_investment_units_milp_constraint(
         self, scenario, exists, override_dict
@@ -1549,7 +1545,7 @@ def test_techs_unit_capacity_max_systemwide_no_transmission_milp_constraint(self
         m.build()
         assert "unit_capacity_max_systemwide_milp" in m.backend.constraints
 
-    @pytest.mark.parametrize("tech", (("test_storage"), ("test_transmission_elec")))
+    @pytest.mark.parametrize("tech", [("test_storage"), ("test_transmission_elec")])
     def test_asynchronous_flow_constraint(self, tech):
         """
         Binary switch for flow in/out can be activated using the option
@@ -1723,7 +1719,7 @@ def simple_supply_longnames(self):
         m.backend.verbose_strings()
         return m
 
-    @pytest.fixture
+    @pytest.fixture()
     def temp_path(self, tmpdir_factory):
         return tmpdir_factory.mktemp("custom_math")
 
@@ -1734,60 +1730,6 @@ def test_new_build_optimal(self, simple_supply):
         assert hasattr(simple_supply, "results")
         assert simple_supply._model_data.attrs["termination_condition"] == "optimal"
 
-    @pytest.mark.parametrize("mode", ["operate", "spores"])
-    def test_add_run_mode_custom_math(self, caplog, mode):
-        caplog.set_level(logging.DEBUG)
-        mode_custom_math = AttrDict.from_yaml(
-            importlib.resources.files("calliope") / "math" / f"{mode}.yaml"
-        )
-        m = build_model({}, "simple_supply,two_hours,investment_costs")
-
-        base_math = deepcopy(m.math)
-        base_math.union(mode_custom_math, allow_override=True)
-
-        backend = PyomoBackendModel(m.inputs, mode=mode)
-        backend._add_run_mode_math()
-
-        assert f"Updating math formulation with {mode} mode math." in caplog.text
-
-        assert m.math != base_math
-        assert backend.inputs.attrs["math"].as_dict() == base_math.as_dict()
-
-    def test_add_run_mode_custom_math_before_build(self, caplog, temp_path):
-        """A user can override the run mode math by including it directly in the additional math list"""
-        caplog.set_level(logging.DEBUG)
-        custom_math = AttrDict({"variables": {"flow_cap": {"active": True}}})
-        file_path = temp_path.join("custom-math.yaml")
-        custom_math.to_yaml(file_path)
-
-        m = build_model(
-            {"config.init.add_math": ["operate", str(file_path)]},
-            "simple_supply,two_hours,investment_costs",
-        )
-        backend = PyomoBackendModel(m.inputs, mode="operate")
-        backend._add_run_mode_math()
-
-        # We set operate mode explicitly in our additional math so it won't be added again
-        assert "Updating math formulation with operate mode math." not in caplog.text
-
-        # operate mode set it to false, then our math set it back to active
-        assert m.math.variables.flow_cap.active
-        # operate mode set it to false and our math did not override that
-        assert not m.math.variables.storage_cap.active
-
-    def test_run_mode_mismatch(self):
-        m = build_model(
-            {"config.init.add_math": ["operate"]},
-            "simple_supply,two_hours,investment_costs",
-        )
-        backend = PyomoBackendModel(m.inputs)
-        with pytest.warns(exceptions.ModelWarning) as excinfo:
-            backend._add_run_mode_math()
-
-        assert check_error_or_warning(
-            excinfo, "Running in plan mode, but run mode(s) {'operate'}"
-        )
-
     @pytest.mark.parametrize(
         "component_type", ["variable", "global_expression", "parameter", "constraint"]
     )
@@ -2131,7 +2073,7 @@ def test_delete_inexistent_pyomo_list(self, simple_supply, component):
         assert "bar" not in getattr(backend_instance, component).keys()
 
     @pytest.mark.parametrize(
-        ["component", "eq"],
+        ("component", "eq"),
         [("global_expressions", "flow_cap + 1"), ("constraints", "flow_cap >= 1")],
     )
     def test_add_allnull_expr_or_constr(self, simple_supply, component, eq):
@@ -2212,7 +2154,7 @@ def test_object_string_representation(self, simple_supply):
         assert not simple_supply.backend.variables.flow_out.coords_in_name
 
     @pytest.mark.parametrize(
-        ["objname", "dims", "objtype"],
+        ("objname", "dims", "objtype"),
         [
             (
                 "flow_out",
@@ -2485,31 +2427,31 @@ def test_has_integer_or_binary_variables_milp(self, supply_milp):
 
 
 class TestShadowPrices:
-    @pytest.fixture(scope="function")
+    @pytest.fixture()
     def simple_supply(self):
         m = build_model({}, "simple_supply,two_hours,investment_costs")
         m.build()
         return m
 
-    @pytest.fixture(scope="function")
+    @pytest.fixture()
     def supply_milp(self):
         m = build_model({}, "supply_milp,two_hours,investment_costs")
         m.build()
         return m
 
-    @pytest.fixture(scope="function")
+    @pytest.fixture()
     def simple_supply_with_yaml_shadow_prices(self):
         m = build_model({}, "simple_supply,two_hours,investment_costs,shadow_prices")
         m.build()
         return m
 
-    @pytest.fixture(scope="function")
+    @pytest.fixture()
     def simple_supply_yaml(self):
         m = build_model({}, "simple_supply,two_hours,investment_costs,shadow_prices")
         m.build()
         return m
 
-    @pytest.fixture(scope="function")
+    @pytest.fixture()
     def simple_supply_yaml_invalid(self):
         m = build_model(
             {},
@@ -2518,7 +2460,7 @@ def simple_supply_yaml_invalid(self):
         m.build()
         return m
 
-    @pytest.fixture(scope="function")
+    @pytest.fixture()
     def supply_milp_yaml(self):
         m = build_model({}, "supply_milp,two_hours,investment_costs,shadow_prices")
         m.build()
diff --git a/tests/test_core_model.py b/tests/test_core_model.py
index 15d34656..67ea9791 100644
--- a/tests/test_core_model.py
+++ b/tests/test_core_model.py
@@ -1,3 +1,4 @@
+import importlib
 import logging
 
 import calliope
@@ -11,6 +12,17 @@
 LOGGER = "calliope.model"
 
 
+@pytest.fixture(scope="module")
+def temp_path(tmpdir_factory):
+    return tmpdir_factory.mktemp("custom_math")
+
+
+@pytest.fixture(scope="module")
+def initialised_model():
+    m = build_model({}, "simple_supply,two_hours,investment_costs")
+    return m
+
+
 class TestModel:
     @pytest.fixture(scope="module")
     def national_scale_example(self):
@@ -30,29 +42,25 @@ def test_info_simple_model(self, simple_supply):
         simple_supply.info()
 
 
-class TestAddMath:
-    @pytest.fixture
-    def storage_inter_cluster(self):
-        return build_model(
-            {"config.init.add_math": ["storage_inter_cluster"]},
-            "simple_supply,two_hours,investment_costs",
+class TestInitMath:
+    def test_apply_base_math_by_default(self, initialised_model):
+        """Base math should be applied by default."""
+        _, applied_math = initialised_model._math_init_from_yaml(
+            initialised_model.config.init
         )
+        assert "base" in applied_math
+        assert not initialised_model.config.init.add_math
 
-    @pytest.fixture
-    def temp_path(self, tmpdir_factory):
-        return tmpdir_factory.mktemp("custom_math")
-
-    def test_internal_override(self, storage_inter_cluster):
-        assert "storage_intra_max" in storage_inter_cluster.math["constraints"].keys()
-
-    def test_variable_bound(self, storage_inter_cluster):
-        assert (
-            storage_inter_cluster.math["variables"]["storage"]["bounds"]["min"]
-            == -np.inf
-        )
+    def test_base_math_deactivation(self, initialised_model):
+        """Base math should be deactivated if requested in config."""
+        init_config = initialised_model.config.init.copy()
+        init_config["base_math"] = False
+        math, applied_math = initialised_model._math_init_from_yaml(init_config)
+        assert not math
+        assert not applied_math
 
     @pytest.mark.parametrize(
-        ["override", "expected"],
+        ("add_math", "in_error"),
         [
             (["foo"], ["foo"]),
             (["bar", "foo"], ["bar", "foo"]),
@@ -60,18 +68,17 @@ def test_variable_bound(self, storage_inter_cluster):
             (["foo.yaml"], ["foo.yaml"]),
         ],
     )
-    def test_allowed_internal_constraint(self, override, expected):
+    def test_math_loading_invalid(self, initialised_model, add_math, in_error):
+        init_config = initialised_model.config.init.copy()
+        init_config.add_math = add_math
         with pytest.raises(calliope.exceptions.ModelError) as excinfo:
-            build_model(
-                {"config.init.add_math": override},
-                "simple_supply,two_hours,investment_costs",
-            )
+            initialised_model._math_init_from_yaml(init_config)
         assert check_error_or_warning(
             excinfo,
-            f"Attempted to load additional math that does not exist: {expected}",
+            f"Attempted to load additional math that does not exist: {in_error}",
         )
 
-    def test_internal_override_from_yaml(self, temp_path):
+    def test_internal_override_from_yaml(self, initialised_model, temp_path):
         new_constraint = calliope.AttrDict(
             {
                 "constraints": {
@@ -83,12 +90,32 @@ def test_internal_override_from_yaml(self, temp_path):
                 }
             }
         )
-        new_constraint.to_yaml(temp_path.join("custom-math.yaml"))
-        m = build_model(
-            {"config.init.add_math": [str(temp_path.join("custom-math.yaml"))]},
+        added_math_path = temp_path.join("custom-math.yaml")
+        new_constraint.to_yaml(added_math_path)
+        init_config = initialised_model.config.init
+        init_config["add_math"] = [str(added_math_path)]
+
+        math, applied_math = initialised_model._math_init_from_yaml(init_config)
+        assert "constraint_name" in math.constraints.keys()
+        assert str(added_math_path) in applied_math
+
+
+class TestModelMathOverrides:
+    @pytest.fixture()
+    def storage_inter_cluster(self):
+        return build_model(
+            {"config.init.add_math": ["storage_inter_cluster"]},
             "simple_supply,two_hours,investment_costs",
         )
-        assert "constraint_name" in m.math["constraints"].keys()
+
+    def test_internal_override(self, storage_inter_cluster):
+        assert "storage_intra_max" in storage_inter_cluster.math["constraints"].keys()
+
+    def test_variable_bound(self, storage_inter_cluster):
+        assert (
+            storage_inter_cluster.math["variables"]["storage"]["bounds"]["min"]
+            == -np.inf
+        )
 
     def test_override_existing_internal_constraint(self, temp_path, simple_supply):
         file_path = temp_path.join("custom-math.yaml")
@@ -180,7 +207,7 @@ def test_base_math(self, caplog, simple_supply):
         ]
 
     @pytest.mark.parametrize(
-        ["equation", "where"],
+        ("equation", "where"),
         [
             ("1 == 1", "True"),
             (
@@ -277,3 +304,37 @@ def test_operate_timeseries(self, operate_model):
             operate_model.results.timesteps
             == pd.date_range("2005-01", "2005-01-02 23:00:00", freq="h")
         )
+
+
+class TestMathUpdate:
+
+    TEST_MODES = ["operate", "spores"]
+
+    @pytest.mark.parametrize("mode", TEST_MODES)
+    def test_math_addition(self, initialised_model, mode):
+        """Run mode math must be added to the model if not present."""
+        mode_custom_math = calliope.AttrDict.from_yaml(
+            importlib.resources.files("calliope") / "math" / f"{mode}.yaml"
+        )
+        new_math = initialised_model.math.copy()
+        new_math.union(mode_custom_math, allow_override=True)
+
+        updated_math, applied_math = initialised_model._math_update_with_mode(
+            {**initialised_model.config.build, **{"mode": mode}}
+        )
+        assert new_math == updated_math
+        assert mode in applied_math
+
+    def test_no_update(self, initialised_model):
+        updated_math, applied_math = initialised_model._math_update_with_mode(
+            initialised_model.config.build
+        )
+        assert updated_math is None
+        assert "base" in applied_math
+
+    @pytest.mark.parametrize("mode", TEST_MODES)
+    def test_mismatch_warning(self, mode):
+        """Warn users is they load unused pre-defined math."""
+        m = build_model({}, "simple_supply,two_hours,investment_costs", add_math=[mode])
+        with pytest.warns(calliope.exceptions.ModelWarning):
+            m._math_update_with_mode(m.config.build)