diff --git a/src/dagmc.cpp b/src/dagmc.cpp
index 56722fde66..6a4865c391 100644
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@@ -280,6 +280,10 @@ void DAGUniverse::init_geometry()
     s->id_ = adjust_geometry_ids_ ? next_surf_id++
                                   : dagmc_instance_->id_by_index(2, i + 1);
 
+    // set surface source attribute if needed
+    if (contains(settings::source_write_surf_id, s->id_))
+      s->surf_source_ = true;
+
     // set BCs
     std::string bc_value =
       dmd_ptr->get_surface_property("boundary", surf_handle);
diff --git a/tests/regression_tests/dagmc/legacy/test.py b/tests/regression_tests/dagmc/legacy/test.py
index 7e4e6e1340..884264f30e 100644
--- a/tests/regression_tests/dagmc/legacy/test.py
+++ b/tests/regression_tests/dagmc/legacy/test.py
@@ -1,9 +1,13 @@
+from pathlib import Path
+
 import openmc
 import openmc.lib
 
-from pathlib import Path
+import h5py
+import numpy as np
 import pytest
-from tests.testing_harness import PyAPITestHarness
+
+from tests.testing_harness import PyAPITestHarness, config
 
 pytestmark = pytest.mark.skipif(
     not openmc.lib._dagmc_enabled(),
@@ -13,7 +17,7 @@ pytestmark = pytest.mark.skipif(
 def model():
     openmc.reset_auto_ids()
 
-    model = openmc.model.Model()
+    model = openmc.Model()
 
     # settings
     model.settings.batches = 5
@@ -73,6 +77,32 @@ def test_missing_material_name(model):
     assert exp_error_msg in str(exec_info.value)
 
 
+def test_surf_source(model):
+    # create a surface source read on this model to ensure
+    # particles are being generated correctly
+    n = 100
+    model.settings.surf_source_write = {'surface_ids': [1], 'max_particles': n}
+
+    # If running in MPI mode, setup proper keyword arguments for run()
+    kwargs = {'openmc_exec': config['exe']}
+    if config['mpi']:
+        kwargs['mpi_args'] = [config['mpiexec'], '-n', config['mpi_np']]
+    model.run(**kwargs)
+
+    with h5py.File('surface_source.h5') as fh:
+        assert fh.attrs['filetype'] == b'source'
+        arr = fh['source_bank'][...]
+    expected_size = n * int(config['mpi_np']) if config['mpi'] else n
+    assert arr.size == expected_size
+
+    # check that all particles are on surface 1 (radius = 7)
+    xs = arr[:]['r']['x']
+    ys = arr[:]['r']['y']
+    rad = np.sqrt(xs**2 + ys**2)
+    assert np.allclose(rad, 7.0)
+
+
 def test_dagmc(model):
     harness = PyAPITestHarness('statepoint.5.h5', model)
     harness.main()
+