diff --git a/CHANGELOG.md b/CHANGELOG.md
index 12c6ad8c84..a9d9b4fa92 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Bug fixes
 
+- Fixed a bug where coordinate systems of variables do not get checked against known ([#3394](https://github.com/pybamm-team/PyBaMM/pull/3394))
 - Fixed a bug where simulations using the CasADi-based solvers would fail randomly with the half-cell model ([#3494](https://github.com/pybamm-team/PyBaMM/pull/3494))
 - Fixed bug that made identical Experiment steps with different end times crash ([#3516](https://github.com/pybamm-team/PyBaMM/pull/3516))
 - Fixed bug in calculation of theoretical energy that made it very slow ([#3506](https://github.com/pybamm-team/PyBaMM/pull/3506))
diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index 2f30da9a5e..4b887a82a7 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -132,6 +132,11 @@ def __init__(
             raise pybamm.DomainError(
                 "domain cannot be particle if name is '{}'".format(name)
             )
+        # check coord_sys against KNOWN_COORD_SYS
+        if coord_sys is not None and coord_sys not in KNOWN_COORD_SYS:
+            raise ValueError(
+                f"Coordinate system is {coord_sys}, not in"
+                " cartesian, cylindrical polar, spherical polar, None")
 
     def create_copy(self):
         """See :meth:`pybamm.Symbol.new_copy()`."""
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index ecbae69796..8834ad716d 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -135,22 +135,30 @@ def divergence(self, symbol, discretised_symbol, boundary_conditions):
         submesh = self.mesh[symbol.domain]
 
         divergence_matrix = self.divergence_matrix(symbol.domains)
+        coord_sys = submesh.coord_sys
 
         # check coordinate system
-        if submesh.coord_sys in ["cylindrical polar", "spherical polar"]:
+        if coord_sys in ["cylindrical polar", "spherical polar"]:
             second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains)
             # create np.array of repeated submesh.edges
             r_edges_numpy = np.kron(np.ones(second_dim_repeats), submesh.edges)
             r_edges = pybamm.Vector(r_edges_numpy)
-            if submesh.coord_sys == "spherical polar":
+            if coord_sys == "spherical polar":
                 out = divergence_matrix @ ((r_edges**2) * discretised_symbol)
-            elif submesh.coord_sys == "cylindrical polar":
+            elif coord_sys == "cylindrical polar":
                 out = divergence_matrix @ (r_edges * discretised_symbol)
-        else:
+        elif coord_sys == "cartesian":
             out = divergence_matrix @ discretised_symbol
-
+        else:
+            if coord_sys is not None:
+                raise ValueError(
+                    f"Coordinate system is {coord_sys}, not in"
+                    " cartesian, cylindrical polar, spherical polar, None"
+                )
         return out
 
+
+
     def divergence_matrix(self, domains):
         """
         Divergence matrix for finite volumes in the appropriate domain.
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 91a5b70044..be7075d566 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
@@ -74,6 +74,36 @@ def test_concatenation(self):
         with self.assertRaisesRegex(pybamm.ShapeError, "child must have size n_nodes"):
             fin_vol.concatenation(edges)
 
+    def test_invalid_coordinate_system(self):
+        with self.assertRaises(ValueError):
+            pybamm.SpatialVariableEdge(
+                "x_s",
+                domain=["separator"],
+                auxiliary_domains={"secondary": "current collector"},
+                coord_sys="cartesian_2",
+            )
+
+    def test_divergence(self):
+        """ unsuccessful try:
+        from unittest.mock import Mock, patch
+        variable = Mock()
+        variable.domain = "some_domain"
+
+        # Mock the necessary parts of the SpatialMethod class
+        spatial_method = pybamm.SpatialMethod()
+        spatial_method.mesh = {"some_domain": Mock()}
+        spatial_method.mesh[
+            "some_domain"].coord_sys = "unsupported_coord_system"
+
+        # Test that the ValueError is raised for unsupported coordinate system
+        with self.assertRaises(ValueError) as context:
+            spatial_method.divergence(variable, None, None)
+
+        self.assertIn(
+            "Coordinate system is unsupported_coord_system", 
+            str(context.exception)) """
+        return 
+
     def test_discretise_diffusivity_times_spatial_operator(self):
         # Setup mesh and discretisation
         mesh = get_mesh_for_testing()