diff --git a/examples/advection_diffusion_multi_material.py b/examples/advection_diffusion_multi_material.py
index a1869335a..b28d53a17 100644
--- a/examples/advection_diffusion_multi_material.py
+++ b/examples/advection_diffusion_multi_material.py
@@ -127,11 +127,6 @@ def velocity_func(x):
 
 # we should be able to automate this
 my_model.interfaces = [F.Interface(6, (bottom_domain, top_domain))]
-my_model.surface_to_volume = {
-    top_surface: top_domain,
-    bottom_surface: bottom_domain,
-    inlet_surf: bottom_domain,
-}
 
 H = F.Species("H", mobile=True)
 
diff --git a/examples/multi_material_1d.py b/examples/multi_material_1d.py
index 2fb2a815c..6fbd9ec31 100644
--- a/examples/multi_material_1d.py
+++ b/examples/multi_material_1d.py
@@ -55,10 +55,6 @@
     left_surface,
     right_surface,
 ]
-my_model.surface_to_volume = {
-    right_surface: right_domain,
-    left_surface: left_domain,
-}
 
 H = F.Species("H", mobile=True)
 trapped_H = F.Species("H_trapped", mobile=False)
diff --git a/examples/multi_material_2d.py b/examples/multi_material_2d.py
index 7a64cbe8c..eb7e49e38 100644
--- a/examples/multi_material_2d.py
+++ b/examples/multi_material_2d.py
@@ -85,10 +85,6 @@ def half(x):
 
 # we should be able to automate this
 my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))]
-my_model.surface_to_volume = {
-    top_surface: top_domain,
-    bottom_surface: bottom_domain,
-}
 
 H = F.Species("H", mobile=True)
 
diff --git a/examples/multi_material_2d_bis.py b/examples/multi_material_2d_bis.py
index cc6ad33da..384052a94 100644
--- a/examples/multi_material_2d_bis.py
+++ b/examples/multi_material_2d_bis.py
@@ -85,10 +85,6 @@ def half(x):
 
 # we should be able to automate this
 my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))]
-my_model.surface_to_volume = {
-    top_surface: top_domain,
-    bottom_surface: bottom_domain,
-}
 
 H = F.Species("H", mobile=True)
 trapped_H = F.Species("H_trapped", mobile=False)
diff --git a/examples/multi_material_transient.py b/examples/multi_material_transient.py
index 8495336fa..a98c5e274 100644
--- a/examples/multi_material_transient.py
+++ b/examples/multi_material_transient.py
@@ -50,10 +50,6 @@
     left_surface,
     right_surface,
 ]
-my_model.surface_to_volume = {
-    right_surface: right_domain,
-    left_surface: left_domain,
-}
 
 H = F.Species("H", mobile=True)
 trapped_H = F.Species("H_trapped", mobile=False)
diff --git a/examples/multi_material_with_one_material.py b/examples/multi_material_with_one_material.py
index ca8c80422..d58f3ffc2 100644
--- a/examples/multi_material_with_one_material.py
+++ b/examples/multi_material_with_one_material.py
@@ -22,7 +22,6 @@
 right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])
 
 my_model.subdomains = [subdomain, left_surface, right_surface]
-my_model.surface_to_volume = {right_surface: subdomain, left_surface: subdomain}
 
 H = F.Species("H", mobile=True)
 trapped_H = F.Species("H_trapped", mobile=False)
diff --git a/src/festim/__init__.py b/src/festim/__init__.py
index 823d519d9..b6133ed86 100644
--- a/src/festim/__init__.py
+++ b/src/festim/__init__.py
@@ -79,5 +79,6 @@
     VolumeSubdomain,
     VolumeSubdomain1D,
     find_volume_from_id,
+    map_surface_to_volume_subdomains,
 )
 from .trap import Trap
diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 9eefe21ab..057683ade 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1055,7 +1055,6 @@ def __init__(
         exports=None,
         traps=None,
         interfaces: list[_subdomain.Interface] | None = None,
-        surface_to_volume: dict | None = None,
         petsc_options: dict | None = None,
     ):
         """Class for a multi-material hydrogen transport problem
@@ -1093,7 +1092,7 @@ def __init__(
             petsc_options=petsc_options,
         )
         self.interfaces = interfaces or []
-        self.surface_to_volume = surface_to_volume or {}
+        self.surface_to_volume = {}
         self.subdomain_to_species = {}  # maps subdomain to species defined in it
 
     @property
@@ -1141,6 +1140,27 @@ def initialise(self):
                 raise TypeError("subdomains attribute in Species should be list")
 
         self.define_meshtags_and_measures()
+        if self.surface_to_volume:
+            # tell users that this is no longer required
+            warnings.warn(
+                f"The surface_to_volume attribute of the {self.__class__.__name__}"
+                " class is no longer required and can be removed."
+                "The mapping between surface and volume subdomains is now done"
+                "automatically based on the connectivity of the mesh and the meshtags",
+                DeprecationWarning,
+            )
+        else:
+            facet_to_cell = self.mesh.mesh.topology.connectivity(
+                self.mesh.mesh.topology.dim - 1, self.mesh.mesh.topology.dim
+            )
+            self.surface_to_volume = _subdomain.map_surface_to_volume_subdomains(
+                ft=self.facet_meshtags,
+                ct=self.volume_meshtags,
+                facet_to_cell=facet_to_cell,
+                volume_subdomains=self.volume_subdomains,
+                surface_subdomains=self.surface_subdomains,
+                comm=self.mesh.mesh.comm,
+            )
 
         # create submeshes and transfer meshtags to subdomains
         for subdomain in self.volume_subdomains:
diff --git a/src/festim/subdomain/__init__.py b/src/festim/subdomain/__init__.py
index d3fb7ec57..7679593f2 100644
--- a/src/festim/subdomain/__init__.py
+++ b/src/festim/subdomain/__init__.py
@@ -1,6 +1,10 @@
 from .interface import Interface
 from .surface_subdomain import SurfaceSubdomain, SurfaceSubdomain1D
-from .volume_subdomain import VolumeSubdomain, VolumeSubdomain1D
+from .volume_subdomain import (
+    VolumeSubdomain,
+    VolumeSubdomain1D,
+    map_surface_to_volume_subdomains,
+)
 
 __all__ = [
     "Interface",
@@ -9,4 +13,5 @@
     "SurfaceSubdomain1D",
     "VolumeSubdomain",
     "VolumeSubdomain1D",
+    "map_surface_to_volume_subdomains",
 ]
diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py
index 22a735b00..5fcc40269 100644
--- a/src/festim/subdomain/volume_subdomain.py
+++ b/src/festim/subdomain/volume_subdomain.py
@@ -9,6 +9,7 @@
 from scifem.mesh import transfer_meshtags_to_submesh
 
 from festim.material import Material
+from festim.subdomain.surface_subdomain import SurfaceSubdomain
 
 # Define the appropriate method based on the version
 try:
@@ -185,3 +186,89 @@ def find_volume_from_id(id: int, volumes: list):
         if vol.id == id:
             return vol
     raise ValueError(f"id {id} not found in list of volumes")
+
+
+def map_surface_to_volume_subdomains(
+    ft: dolfinx.mesh.MeshTags,
+    ct: dolfinx.mesh.MeshTags,
+    facet_to_cell: dolfinx.cpp.graph.AdjacencyList_int32,
+    volume_subdomains: list[VolumeSubdomain],
+    surface_subdomains: list[SurfaceSubdomain],
+    comm=None,
+) -> dict[SurfaceSubdomain, VolumeSubdomain]:
+    """Maps surface subdomains to volume subdomains based on the facet and cell meshtags
+    and the facet to cell connectivity.
+
+
+    Raises:
+        AssertionError: if a surface subdomain is connected to multiple volume subdomains
+
+    Args:
+        ft: the facet meshtags of the parent mesh
+        ct: the cell meshtags of the parent mesh
+        facet_to_cell: the facet to cell connectivity of the parent mesh
+        volume_subdomains: the list of volume subdomains
+        surface_subdomains: the list of surface subdomains
+        comm: MPI communicator (required for parallel runs)
+
+    Returns:
+        dict[SurfaceSubdomain, VolumeSubdomain]: a dictionary mapping surface subdomains
+            to volume subdomains
+    """
+
+    # get connected cells for tagged facets
+    start_indices = facet_to_cell.offsets[ft.indices]
+    end_indices = facet_to_cell.offsets[ft.indices + 1]
+    num_connections = end_indices - start_indices
+
+    # A facet is connected to at most 2 cells (boundary = 1, interior = 2)
+    cell_ids_0 = facet_to_cell.array[start_indices]
+    has_second_cell = num_connections == 2
+    cell_ids_1 = facet_to_cell.array[start_indices[has_second_cell] + 1]
+
+    connected_cells = np.concatenate([cell_ids_0, cell_ids_1])
+    connected_facet_tags = np.concatenate([ft.values, ft.values[has_second_cell]])
+
+    # map connected cells to their cell tags
+    sort_idx = np.argsort(ct.indices)
+    sorted_ct_indices = ct.indices[sort_idx]
+    sorted_ct_values = ct.values[sort_idx]
+
+    idx = np.searchsorted(sorted_ct_indices, connected_cells)
+    # mask out-of-bounds
+    valid = idx < len(sorted_ct_indices)
+    # of those in bounds, check if they actually match
+    valid[valid] = sorted_ct_indices[idx[valid]] == connected_cells[valid]
+
+    valid_cell_tags = sorted_ct_values[idx[valid]]
+    valid_facet_tags = connected_facet_tags[valid]
+
+    unique_pairs = np.unique(np.vstack((valid_facet_tags, valid_cell_tags)).T, axis=0)
+    if comm is not None and comm.size > 1:
+        all_pairs = comm.allgather(unique_pairs)
+        non_empty = [p for p in all_pairs if len(p) > 0]
+        if non_empty:
+            unique_pairs = np.unique(np.vstack(non_empty), axis=0)
+
+    surface_tag_to_subdomain = {s.id: s for s in surface_subdomains}
+    volume_tag_to_subdomain = {v.id: v for v in volume_subdomains}
+
+    surface_to_subdomain = {}
+
+    for s_tag, v_tag in unique_pairs:
+        dolfinx.log.log(
+            dolfinx.log.LogLevel.INFO,
+            f"Facet tag {s_tag} is connected to cell tag {v_tag}",
+        )
+        s_subdomain = surface_tag_to_subdomain.get(s_tag)
+        v_subdomain = volume_tag_to_subdomain.get(v_tag)
+
+        if s_subdomain and v_subdomain:
+            if s_subdomain in surface_to_subdomain:
+                assert surface_to_subdomain[s_subdomain] == v_subdomain, (
+                    f"Surface subdomain {s_subdomain.id} is connected to multiple volume subdomains: "
+                    f"{surface_to_subdomain[s_subdomain].id} and {v_subdomain.id}"
+                )
+            else:
+                surface_to_subdomain[s_subdomain] = v_subdomain
+    return surface_to_subdomain
diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py
index 3d3e728a0..3a11b1fea 100644
--- a/test/system_tests/test_advection_term_integration.py
+++ b/test/system_tests/test_advection_term_integration.py
@@ -250,11 +250,6 @@ def velocity_func(x):
 
     # we should be able to automate this
     my_model.interfaces = [F.Interface(6, (bottom_domain, top_domain))]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-        inlet_surf: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
diff --git a/test/system_tests/test_cylindrical_coordinates.py b/test/system_tests/test_cylindrical_coordinates.py
index eaa313e98..460fa1249 100644
--- a/test/system_tests/test_cylindrical_coordinates.py
+++ b/test/system_tests/test_cylindrical_coordinates.py
@@ -115,10 +115,6 @@ def c_exact_right(x):
     my_model.subdomains = [vol_1, vol_2, left, right]
 
     my_model.interfaces = [F.Interface(5, (vol_1, vol_2), penalty_term=100)]
-    my_model.surface_to_volume = {
-        left: vol_1,
-        right: vol_2,
-    }
 
     H = F.Species("H", mobile=True, subdomains=[vol_1, vol_2])
     my_model.species = [H]
diff --git a/test/system_tests/test_misc.py b/test/system_tests/test_misc.py
index 47d855942..5b9152017 100644
--- a/test/system_tests/test_misc.py
+++ b/test/system_tests/test_misc.py
@@ -210,11 +210,6 @@ def test_temp_dependent_bc_mixed_domain_temperature_as_function():
         bottom_surface,
     ]
 
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
-
     H = F.Species("H", mobile=True, subdomains=[bottom_domain, top_domain])
 
     my_model.species = [H]
diff --git a/test/system_tests/test_multi_mat_penalty.py b/test/system_tests/test_multi_mat_penalty.py
index 116da435e..f5ef5c132 100644
--- a/test/system_tests/test_multi_mat_penalty.py
+++ b/test/system_tests/test_multi_mat_penalty.py
@@ -112,10 +112,6 @@ def c_exact_bot_np(x):
     my_model.interfaces = [
         F.Interface(5, (bottom_domain, top_domain), penalty_term=1),
     ]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
@@ -195,10 +191,6 @@ def test_derived_quantities_multi_mat():
     my_model.interfaces = [
         F.Interface(5, (bottom_domain, top_domain), penalty_term=1),
     ]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
@@ -270,11 +262,6 @@ def test_penalty_multispecies():
     for spe in my_model.species:
         spe.subdomains = [vol1, vol2]
 
-    my_model.surface_to_volume = {
-        left_surf: vol1,
-        right_surf: vol2,
-    }
-
     my_model.boundary_conditions = [
         # Protium BCs
         F.FixedConcentrationBC(left_surf, value=10, species=protium),
diff --git a/test/system_tests/test_multi_material.py b/test/system_tests/test_multi_material.py
index 3a7c1551a..41a31ee19 100644
--- a/test/system_tests/test_multi_material.py
+++ b/test/system_tests/test_multi_material.py
@@ -110,10 +110,6 @@ def c_exact_bot_np(x):
     ]
 
     my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
@@ -182,10 +178,6 @@ def test_1_material_discontinuous_version(tmpdir):
     right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])
 
     my_model.subdomains = [subdomain, left_surface, right_surface]
-    my_model.surface_to_volume = {
-        right_surface: subdomain,
-        left_surface: subdomain,
-    }
 
     H = F.Species("H", mobile=True)
     trapped_H = F.Species("H_trapped", mobile=False)
@@ -278,10 +270,6 @@ def test_3_materials_transient(tmpdir):
         left_surface,
         right_surface,
     ]
-    my_model.surface_to_volume = {
-        right_surface: right_domain,
-        left_surface: left_domain,
-    }
 
     H = F.Species("H", mobile=True)
     trapped_H = F.Species("H_trapped", mobile=False)
@@ -367,10 +355,6 @@ def test_2_mats_particle_flux_bc(tmpdir):
 
     # we should be able to automate this
     my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
@@ -418,10 +402,6 @@ def test_all_cells_are_not_tagged():
 
     my_model.subdomains = [vol1, vol2, surf1, surf2]
 
-    my_model.surface_to_volume = {
-        surf1: vol1,
-        surf2: vol2,
-    }
     my_model.interfaces = [F.Interface(id=3, subdomains=[vol1, vol2])]
 
     my_model.mesh = F.Mesh1D(np.linspace(0, 1, 10))
diff --git a/test/system_tests/test_spherical_coordinates.py b/test/system_tests/test_spherical_coordinates.py
index d477f484d..2b2887405 100644
--- a/test/system_tests/test_spherical_coordinates.py
+++ b/test/system_tests/test_spherical_coordinates.py
@@ -116,10 +116,6 @@ def c_exact_right(x):
     my_model.subdomains = [vol_1, vol_2, left, right]
 
     my_model.interfaces = [F.Interface(5, (vol_1, vol_2), penalty_term=1e5)]
-    my_model.surface_to_volume = {
-        left: vol_1,
-        right: vol_2,
-    }
 
     H = F.Species("H", mobile=True, subdomains=[vol_1, vol_2])
     my_model.species = [H]
diff --git a/test/test_complex_reaction.py b/test/test_complex_reaction.py
index bcbbb7923..e6fbabf18 100644
--- a/test/test_complex_reaction.py
+++ b/test/test_complex_reaction.py
@@ -188,10 +188,6 @@ def test_surface_reaction_in_discontinuous_case():
     right_surface = F.SurfaceSubdomain1D(id=2, x=1)
 
     my_model.subdomains = [subdomain, left_surface, right_surface]
-    my_model.surface_to_volume = {
-        right_surface: subdomain,
-        left_surface: subdomain,
-    }
 
     H = F.Species("H", mobile=True, subdomains=[subdomain])
     my_model.species = [H]
diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py
index d33e002a2..bc41336de 100644
--- a/test/test_h_transport_problem.py
+++ b/test/test_h_transport_problem.py
@@ -1411,11 +1411,6 @@ def test_surface_reaction_BC_discontinuous():
         outlet,
     ]
 
-    my_model.surface_to_volume = {
-        inlet: left_vol,
-        outlet: right_vol,
-    }
-
     my_model.method_interface = F.InterfaceMethod.nitsche
     my_model.interfaces = [
         F.Interface(id=5, subdomains=[left_vol, right_vol], penalty_term=10),
diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py
index 588785445..6214fe278 100644
--- a/test/test_initial_condition.py
+++ b/test/test_initial_condition.py
@@ -347,10 +347,6 @@ def test_initial_condition_discontinuous():
     my_model.interfaces = [
         F.Interface(id=3, subdomains=[vol1, vol2], penalty_term=1000)
     ]
-    my_model.surface_to_volume = {
-        left_surf: vol1,
-        right_surf: vol2,
-    }
 
     my_model.temperature = 300
 
diff --git a/test/test_profile.py b/test/test_profile.py
index c52de4544..bd36291fd 100644
--- a/test/test_profile.py
+++ b/test/test_profile.py
@@ -144,11 +144,6 @@ def test_profile_discontinuous():
     for spe in my_model.species:
         spe.subdomains = [vol1, vol2]
 
-    my_model.surface_to_volume = {
-        left_surf: vol1,
-        right_surf: vol2,
-    }
-
     my_model.boundary_conditions = [
         # Protium BCs
         F.FixedConcentrationBC(left_surf, value=10, species=protium),
@@ -215,11 +210,6 @@ def test_profile_discontinuous_single_species():
     for spe in my_model.species:
         spe.subdomains = [vol1, vol2]
 
-    my_model.surface_to_volume = {
-        left_surf: vol1,
-        right_surf: vol2,
-    }
-
     my_model.boundary_conditions = [
         F.FixedConcentrationBC(left_surf, value=10, species=protium),
         F.FixedConcentrationBC(right_surf, value=0, species=protium),
diff --git a/test/test_subdomains.py b/test/test_subdomains.py
index 351ac0e42..844b0fd38 100644
--- a/test/test_subdomains.py
+++ b/test/test_subdomains.py
@@ -234,3 +234,49 @@ def test_all_cells_are_not_tagged():
         AssertionError, match="All cells must be tagged with a non-zero value"
     ):
         mesh.define_meshtags(surface_subdomains=[], volume_subdomains=[vol1, vol2])
+
+
+def test_map_surface_to_volume_subdomains():
+    """
+    Tests that the function map_surface_to_volume_subdomains correctly maps surface
+    subdomains to volume subdomains based on the facet to cell connectivity
+    """
+
+    n = 20
+    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n)
+
+    festim_mesh = F.Mesh(mesh)
+
+    surface_1 = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 0))
+    surface_2 = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1))
+
+    material_1 = F.Material(D_0=1, E_D=0, name="material_1")
+    volume_1 = F.VolumeSubdomain(
+        id=2, material=material_1, locator=lambda x: x[0] <= 0.5
+    )
+    volume_2 = F.VolumeSubdomain(
+        id=3, material=material_1, locator=lambda x: x[0] >= 0.5
+    )
+
+    ft, ct = festim_mesh.define_meshtags(
+        surface_subdomains=[surface_1, surface_2],
+        volume_subdomains=[volume_1, volume_2],
+    )
+
+    facet_to_cell = mesh.topology.connectivity(mesh.topology.dim - 1, mesh.topology.dim)
+
+    surface_to_subdomain = F.map_surface_to_volume_subdomains(
+        ft, ct, facet_to_cell, [volume_1, volume_2], [surface_1, surface_2]
+    )
+    print(surface_to_subdomain)
+    for surface, volume in surface_to_subdomain.items():
+        print(f"Surface {surface.id} is connected to volume {volume.id}")
+
+    assert surface_to_subdomain[surface_1] == volume_1, (
+        "Expected surface 1 to be connected to volume 1, "
+        f"got {surface_to_subdomain[surface_1].id}"
+    )
+    assert surface_to_subdomain[surface_2] == volume_2, (
+        "Expected surface 2 to be connected to volume 2, "
+        f"got {surface_to_subdomain[surface_2].id}"
+    )
diff --git a/test/test_vtx.py b/test/test_vtx.py
index edbb2e531..33ce7bf69 100644
--- a/test/test_vtx.py
+++ b/test/test_vtx.py
@@ -68,7 +68,6 @@ def test_vtx_DG(tmpdir):
 
     my_model.temperature = 55
     my_model.subdomains = [s0, s1, l0, l1]
-    my_model.surface_to_volume = {l0: s0, l1: s1}
     # NOTE: Ask Remi why `H` has to live in both s0 and s1
     my_model.species = [
         F.Species("H", subdomains=[s0, s1]),