diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index a4cf0f4978..d7c2018a9e 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -669,9 +669,9 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_MIGRATION:
-      score =
-        1.0 / 6.0 * p.wgt() *
-        ((p.r() - p.r_born()).norm_sq() - (p.r_last() - p.r_born()).norm_sq());
+      score = 1.0 / 6.0 * p.wgt() *
+              (p.r() - p.r_last())
+                .dot((p.r() - p.r_born()) + (p.r_last() - p.r_born()));
       break;
 
     case SCORE_NU_FISSION:
@@ -1772,9 +1772,9 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_MIGRATION:
-      score =
-        1.0 / 6.0 * p.wgt() *
-        ((p.r() - p.r_born()).norm_sq() - (p.r_last() - p.r_born()).norm_sq());
+      score = 1.0 / 6.0 * p.wgt() *
+              (p.r() - p.r_last())
+                .dot((p.r() - p.r_born()) + (p.r_last() - p.r_born()));
       break;
 
     case SCORE_NU_SCATTER:
@@ -2646,7 +2646,7 @@ void score_surface_tally(Particle& p, const vector<int>& tallies)
       // for a further scoring function.
       double score = current * filter_weight;
       for (auto score_index = 0; score_index < tally.scores_.size();
-           ++score_index) {
+        ++score_index) {
 #pragma omp atomic
         tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
       }
@@ -2721,7 +2721,7 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
 
             // Loop over scores.
             for (auto score_index = 0; score_index < tally.scores_.size();
-                 ++score_index) {
+              ++score_index) {
 #pragma omp atomic
               tally.results_(filter_index, score_index, TallyResult::VALUE) +=
                 filter_weight;
diff --git a/tests/unit_tests/test_migration.py b/tests/unit_tests/test_migration.py
new file mode 100644
index 0000000000..a0aa220361
--- /dev/null
+++ b/tests/unit_tests/test_migration.py
@@ -0,0 +1,53 @@
+import openmc
+import pytest
+import numpy as np
+
+from tests import cdtemp
+
+def sphere_model(radius, boundary_type):
+    
+    model = openmc.Model()
+
+    # Material
+    material = openmc.Material(name="Hydrogen")
+    material.add_nuclide("H1", 1.0)
+    material.set_density('g/cm3', 1.0)
+    material.add_s_alpha_beta('c_H_in_H2O')
+
+    # Geometry
+    sphere = openmc.Sphere(r=radius, boundary_type=boundary_type)
+    cell = openmc.Cell(region=-sphere, fill=material)
+    model.geometry = openmc.Geometry([cell])
+
+    # Settings
+    model.settings.particles = 1000
+    model.settings.batches = 5
+    model.settings.run_mode = 'fixed source'
+    model.settings.source = openmc.IndependentSource()
+
+
+    # Tally
+    tally = openmc.Tally()
+    tally.scores = ["migration-area"]
+    model.tallies = [tally]
+
+    return model, tally
+
+def test_reflection_is_equivalent_to_large_model(run_in_tmpdir):
+    openmc.reset_auto_ids()
+    refl_model, tally1 = sphere_model(2.5, "reflective")
+    large_model, tally2 = sphere_model(100, "vacuum")
+
+    with cdtemp():
+        refl_model.run(apply_tally_results=True)
+        mean1, std1 = (tally1.mean.squeeze(),
+                       tally1.std_dev.squeeze())
+    
+    with cdtemp():
+        large_model.run(apply_tally_results=True)
+        mean2, std2 = (tally2.mean.squeeze(),
+                       tally2.std_dev.squeeze())
+    std = (std1**2+std2**2)**0.5
+    diff = np.abs(mean2-mean1)
+        
+    assert diff <= 3*std