diff --git a/src/particle.cpp b/src/particle.cpp
index 10e8f213dd..3dfcd233ef 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -182,10 +182,21 @@ void Particle::from_source(const SourceSite* src)
   parent_nuclide() = src->parent_nuclide;
   delayed_group() = src->delayed_group;
 
-  // Convert signed surface ID to signed index
+  // If the source has a surface ID assigned, attempt to determine the signed
+  // surface index by comparing the particle's direction with the surface
+  // normal. This only works if the particle is on the surface.
   if (src->surf_id != SURFACE_NONE) {
-    int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
-    surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
+    auto it = model::surface_map.find(std::abs(src->surf_id));
+    if (it != model::surface_map.end()) {
+      int i_surface = it->second;
+      const auto& surf = *model::surfaces[i_surface];
+      if (surf.geom_type() == GeometryType::CSG &&
+          std::abs(surf.evaluate(r())) < FP_COINCIDENT) {
+        int index_plus_one = i_surface + 1;
+        surface() =
+          (u().dot(surf.normal(r())) > 0.0) ? index_plus_one : -index_plus_one;
+      }
+    }
   }
 
   wgt_born() = src->wgt_born;
diff --git a/tests/unit_tests/test_surface_source_write.py b/tests/unit_tests/test_surface_source_write.py
index 6f18d32b71..e229fbb72e 100644
--- a/tests/unit_tests/test_surface_source_write.py
+++ b/tests/unit_tests/test_surface_source_write.py
@@ -193,6 +193,40 @@ def test_particle_direction(parameter, run_in_tmpdir, model):
                 assert False
 
 
+def test_surface_source_id(run_in_tmpdir):
+    """Test handling of surface IDs for source particles."""
+
+    # Surface 20 is the surface matching the source in this model. Surface 11 is
+    # absent and surface 10 is reused for an unrelated boundary that does not
+    # contain the source position.
+    sources = []
+    for direction in ((1.0, 0.0, 0.0), (-1.0, 0.0, 0.0)):
+        for source_surface_id in (20, 11, 10):
+            filename = f"surface_source_{len(sources)}.h5"
+            site = openmc.SourceParticle(
+                r=(0.0, 0.0, 0.0), u=direction, surf_id=source_surface_id
+            )
+            openmc.write_source_file([site], filename)
+            sources.append(openmc.FileSource(filename))
+
+    left = openmc.XPlane(-1.0, surface_id=1, boundary_type="vacuum")
+    middle = openmc.XPlane(0.0, surface_id=20)
+    right = openmc.XPlane(1.0, surface_id=30, boundary_type="vacuum")
+    bottom = openmc.YPlane(-1.0, surface_id=40, boundary_type="vacuum")
+    top = openmc.YPlane(1.0, surface_id=10, boundary_type="vacuum")
+    cells = [
+        openmc.Cell(region=+left & -middle & +bottom & -top),
+        openmc.Cell(region=+middle & -right & +bottom & -top),
+    ]
+
+    model = openmc.Model(geometry=openmc.Geometry(cells))
+    model.settings.run_mode = "fixed source"
+    model.settings.particles = 1000
+    model.settings.batches = 1
+    model.settings.source = sources
+    model.run()
+
+
 @pytest.fixture
 def model_dagmc(request):
     """Model based on the mesh file 'dagmc.h5m' available from