diff --git a/avaframe/com1DFA/DFAToolsCython.pyx b/avaframe/com1DFA/DFAToolsCython.pyx
index 63708e754..809e12d1c 100644
--- a/avaframe/com1DFA/DFAToolsCython.pyx
+++ b/avaframe/com1DFA/DFAToolsCython.pyx
@@ -573,8 +573,8 @@ cpdef (double, double, int, int, int, double, double, double, double) samosProje
     # are we in the same cell?
     if Lxi==Lxj and Lyi==Lyj:
       return xNew, yNew, iCell, Lx0, Ly0, w[0], w[1], w[2], w[3]
-    Lxj = Lxi
-    Lyj = Lyi
+    Lxi = Lxj
+    Lyi = Lyj
   return xNew, yNew, iCell, Lx0, Ly0, w[0], w[1], w[2], w[3]
 
 
@@ -675,6 +675,9 @@ cpdef (double, double, double, int, int, int, double, double, double, double) di
     # second step: conserve the distance dist (correction ov the position)
     # measure distance between this new point and the previous time step position
     distn = norm(xNew-xPrev, yNew-yPrev, zNew-zPrev)
+    # guard against stopped particle (distn == 0)
+    if distn == 0:
+      break
     # adjust position on the Xprev, Xnew line
     xNew = xPrev + (xNew-xPrev) * dist / distn
     yNew = yPrev + (yNew-yPrev) * dist / distn
diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx
index 919707a3c..568eb1268 100644
--- a/avaframe/com1DFA/DFAfunctionsCython.pyx
+++ b/avaframe/com1DFA/DFAfunctionsCython.pyx
@@ -481,7 +481,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
     if entrMassCell < 0:
       entrMassCell = 0
     entrMassRaster[indCellY, indCellX] = entrMassCell
-  fields['entrMassRaster'] = np.asarray(entrMassRaster)
+  fields['entrMassRaster'] = np.asarray(entrMassRaster).copy()
 
   return particles, force, fields
 
@@ -1115,7 +1115,6 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0):
 
   if stop:
     particles['iterate'] = False
-    particles['stoppedParticles'] = particles
     massStopped = massStopped + massFlowing
 
     if typeStop == 1:
@@ -1398,27 +1397,27 @@ def updateFieldsC(cfg, particles, dem, fields):
       FTStopBilinear[j, i] = - MassStopBilinear[j, i] / (areaRaster[j, i] * rho)  # / m * FTBilinear[j, i] 
       FTEntBilinear[j, i] = - MassEntBilinear[j, i] / (areaRaster[j, i] * rhoEnt)
 
-  fields['FM'] = np.asarray(MassBilinear)
-  fields['FV'] = np.asarray(VBilinear)
-  fields['Vx'] = np.asarray(VXBilinear)
-  fields['Vy'] = np.asarray(VYBilinear)
-  fields['Vz'] = np.asarray(VZBilinear)
-  fields['FT'] = np.asarray(FTBilinear)
-  fields['pfv'] = np.asarray(PFV)
-  fields['pft'] = np.asarray(PFT)
-  fields['dmDet'] = np.asarray(DMDet)
-  fields['FTStop'] = np.asarray(FTStopBilinear)
-  fields['FTDet'] = np.asarray(FTDetBilinear)
-  fields['FTEnt'] = np.asarray(FTEntBilinear)
+  fields['FM'] = np.asarray(MassBilinear).copy()
+  fields['FV'] = np.asarray(VBilinear).copy()
+  fields['Vx'] = np.asarray(VXBilinear).copy()
+  fields['Vy'] = np.asarray(VYBilinear).copy()
+  fields['Vz'] = np.asarray(VZBilinear).copy()
+  fields['FT'] = np.asarray(FTBilinear).copy()
+  fields['pfv'] = np.asarray(PFV).copy()
+  fields['pft'] = np.asarray(PFT).copy()
+  fields['dmDet'] = np.asarray(DMDet).copy()
+  fields['FTStop'] = np.asarray(FTStopBilinear).copy()
+  fields['FTDet'] = np.asarray(FTDetBilinear).copy()
+  fields['FTEnt'] = np.asarray(FTEntBilinear).copy()
 
   if computeP:
-    fields['ppr'] = np.asarray(PP)
-    fields['P'] = np.asarray(PBilinear)
+    fields['ppr'] = np.asarray(PP).copy()
+    fields['P'] = np.asarray(PBilinear).copy()
   if computeTA:
-    fields['TA'] = np.asarray(travelAngleField)
-    fields['pta'] = np.asarray(PTA)
+    fields['TA'] = np.asarray(travelAngleField).copy()
+    fields['pta'] = np.asarray(PTA).copy()
   if computeKE:
-    fields['pke'] = np.asarray(PKE)
+    fields['pke'] = np.asarray(PKE).copy()
 
 
   for k in range(nPart):
@@ -1674,7 +1673,9 @@ def computeCohesionForceC(cfg, particles, force):
   force['forceSPHX'] = np.asarray(forceSPHX)
   force['forceSPHY'] = np.asarray(forceSPHY)
   force['forceSPHZ'] = np.asarray(forceSPHZ)
-  particles['bondDist'] = np.asarray(bondDist)
+  # bondDist is modified in-place via the memoryview; do not reassign
+  # particles['bondDist'] here, otherwise numpy will create a VIEW with
+  # base=memoryview which can hang on deallocation on Windows.
   return force, particles
 
 
@@ -1775,9 +1776,9 @@ def initializeBondsC(particles, triangles):
     bondDist[bondStart2[p2+1]-1] = distanceIni
     bondStart2[p2+1] = bondStart2[p2+1] - 1
 
-  particles['bondStart'] = np.asarray(bondStart)
-  particles['bondPart'] = np.asarray(bondPart)
-  particles['bondDist'] = np.asarray(bondDist)
+  particles['bondStart'] = np.asarray(bondStart).copy()
+  particles['bondPart'] = np.asarray(bondPart).copy()
+  particles['bondDist'] = np.asarray(bondDist).copy()
   return particles
 
 def countRemovedBonds(particles, mask, nRemove):
@@ -1870,9 +1871,35 @@ def removedBonds(particles, mask, nRemove, nBondRemove):
          countBond = countBond + 1
 
 
-  particles['bondStart'] = np.asarray(bondStartNew)
-  particles['bondPart'] = np.asarray(bondPartNew)
-  particles['bondDist'] = np.asarray(bondDistNew)
+  # Remap bondStart and bondPart to use compacted (new) particle indices
+  # After removePart compacts particle arrays via np.extract, old indices become stale
+  cdef int[:] oldToNew = np.full(nPart, -1, dtype=np.int32)
+  cdef int newIdx = 0
+  cdef int kk
+  for kk in range(nPart):
+    if keepParticle[kk] == 1:
+      oldToNew[kk] = newIdx
+      newIdx = newIdx + 1
+  cdef int nPartNew = newIdx
+
+  # Rebuild bondStart indexed by new particle indices
+  cdef int[:] bondStartFinal = np.zeros(nPartNew + 1, dtype=np.int32)
+  cdef int bondCumul = 0
+  newIdx = 0
+  for kk in range(nPart):
+    if keepParticle[kk] == 1:
+      bondCumul = bondCumul + (bondStartNew[kk + 1] - bondStartNew[kk])
+      bondStartFinal[newIdx + 1] = bondCumul
+      newIdx = newIdx + 1
+
+  # Remap bondPartNew target indices from old to new
+  cdef int ibb
+  for ibb in range(countBondNew):
+    bondPartNew[ibb] = oldToNew[bondPartNew[ibb]]
+
+  particles['bondStart'] = np.asarray(bondStartFinal).copy()
+  particles['bondPart'] = np.asarray(bondPartNew).copy()
+  particles['bondDist'] = np.asarray(bondDistNew).copy()
   return particles
 
 
diff --git a/avaframe/com1DFA/particleTools.py b/avaframe/com1DFA/particleTools.py
index a8b13fbed..6676407f0 100644
--- a/avaframe/com1DFA/particleTools.py
+++ b/avaframe/com1DFA/particleTools.py
@@ -229,7 +229,8 @@ def removePart(particles, mask, nRemove, reasonString="", snowSlide=0):
     for key in particles:
         if key == "nPart":
             particles["nPart"] = particles["nPart"] - nRemove
-        elif key == "stoppedParticles":
+        elif key == "stoppedParticles" or key == "bondStart" or key == "bondPart" or key == "bondDist":
+            # bond arrays are already rebuilt by removedBonds with new-indexed data
             continue
         # for all keys in particles that are arrays of size nPart do:
         elif type(particles[key]).__module__ == np.__name__:
diff --git a/avaframe/tests/test_DFAToolsCython.py b/avaframe/tests/test_DFAToolsCython.py
index cbcdf0c34..54c4cd1d4 100644
--- a/avaframe/tests/test_DFAToolsCython.py
+++ b/avaframe/tests/test_DFAToolsCython.py
@@ -310,6 +310,105 @@ def test_reprojectionC(capfd):
     assert abs(dist-dist0) <= threshold*(dist0 + csz)
 
 
+def test_samosProjectionConvergence(capfd):
+    """Test samosProjectionIteratrive convergence: result is stable across iteration counts."""
+    ncols = 100
+    nrows = 50
+    csz = 5
+    header = {}
+    header['ncols'] = ncols
+    header['nrows'] = nrows
+    header['cellsize'] = csz
+    dem = {}
+    dem['header'] = header
+    x0 = 300
+    z0 = 500
+    X = np.linspace(0, csz * (ncols - 1), ncols)
+    Y = np.linspace(0, csz * (nrows - 1), nrows)
+    XX, YY = np.meshgrid(X, Y)
+    ZZ = z0 / (x0 * x0) * (XX - x0) * (XX - x0)
+    dem['rasterData'] = ZZ
+    num = 1
+    interpOption = 2
+    dem = gT.getNormalMesh(dem, num=num)
+    Nx = dem['Nx']
+    Ny = dem['Ny']
+    Nz = dem['Nz']
+
+    # Point above the parabola surface, offset to trigger cell change during reprojection
+    xIn = 350
+    yIn = 125
+    Lx0, Ly0, iCell, w0, w1, w2, w3 = DFAfunC.getCellAndWeights(
+        xIn, yIn, ncols, nrows, csz, interpOption)
+    zSurf = DFAfunC.getScalar(Lx0, Ly0, w0, w1, w2, w3, ZZ)
+    zIn = zSurf + 20
+
+    # Reference: many iterations
+    xRef, yRef, iCellRef, LxR0, LyR0, wR0, wR1, wR2, wR3 = DFAfunC.samosProjectionIteratrive(
+        xIn, yIn, zIn, ZZ, Nx, Ny, Nz, csz, ncols, nrows, interpOption, 50)
+    zRef = DFAfunC.getScalar(LxR0, LyR0, wR0, wR1, wR2, wR3, ZZ)
+
+    # With enough iterations, result matches the reference (non-regression)
+    for nIter in [29, 30, 50]:
+        xN, yN, iCellN, LxN0, LyN0, wN0, wN1, wN2, wN3 = DFAfunC.samosProjectionIteratrive(
+            xIn, yIn, zIn, ZZ, Nx, Ny, Nz, csz, ncols, nrows, interpOption, nIter)
+        zN = DFAfunC.getScalar(LxN0, LyN0, wN0, wN1, wN2, wN3, ZZ)
+        assert xN == pytest.approx(xRef, rel=1e-4)
+        assert yN == pytest.approx(yRef, rel=1e-4)
+        assert zN == pytest.approx(zRef, rel=1e-4)
+
+
+def test_distConservStoppedParticle(capfd):
+    """distConservProjectionIteratrive handles distn=0 inside loop without NaN.
+
+    prev is on a flat surface, current is above the surface at the same (x,y).
+    Initial distn is non-zero (current is above surface), so the loop enters.
+    After orthogonal reprojection, the particle lands on the surface and distn
+    becomes zero, triggering the division by zero in the distance conservation step."""
+    ncols = 20
+    nrows = 20
+    csz = 5
+    header = {}
+    header['ncols'] = ncols
+    header['nrows'] = nrows
+    header['cellsize'] = csz
+    dem = {}
+    dem['header'] = header
+    X = np.linspace(0, csz * (ncols - 1), ncols)
+    Y = np.linspace(0, csz * (nrows - 1), nrows)
+    XX, YY = np.meshgrid(X, Y)
+    # Flat surface so normal is straight up, no lateral shift
+    ZZ = np.zeros_like(XX)
+    dem['rasterData'] = ZZ
+    num = 1
+    interpOption = 2
+    dem = gT.getNormalMesh(dem, num=num)
+    Nx = dem['Nx']
+    Ny = dem['Ny']
+    Nz = dem['Nz']
+
+    xPrev = 50.0
+    yPrev = 50.0
+    zPrev = 0.0  # on flat surface
+    xCur = xPrev  # same x, y — no lateral movement
+    yCur = yPrev
+    zCur = 20.0  # above surface
+    threshold = 0.001
+
+    xNew, yNew, zNew, iCell, Lx0, Ly0, w0, w1, w2, w3 = DFAfunC.distConservProjectionIteratrive(
+        xPrev, yPrev, zPrev, ZZ, Nx, Ny, Nz,
+        xCur, yCur, zCur,
+        csz, ncols, nrows, interpOption, 10, threshold)
+
+    # Must not be NaN or inf (would be if division by zero is unguarded)
+    assert not np.isnan(xNew)
+    assert not np.isnan(yNew)
+    assert not np.isnan(zNew)
+    assert np.isfinite(xNew)
+    assert np.isfinite(yNew)
+    assert np.isfinite(zNew)
+
+
 def test_SamosATfric(capfd):
     """ Test the account4FrictionForce function"""
     tau0 = 0
diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py
index b7375191d..123ede57c 100644
--- a/avaframe/tests/test_DFAfunctionsCython.py
+++ b/avaframe/tests/test_DFAfunctionsCython.py
@@ -582,6 +582,112 @@ def test_updatePositionC():
     assert particles["massStopped"] == -np.sum(particles["stoppedParticles"]["m"])
 
 
+def test_updatePositionC_no_circular_reference():
+    """updatePositionC must not create circular reference in stoppedParticles.
+
+    Uses stopCrit=1.0 with typeStop=0 and stopCritType=kinEnergy so that the stop
+    criterion is met on the first call. This exercises the `if stop:` code path
+    where the circular reference was previously created."""
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {
+        "stopCrit": "1.01",
+        "stopCritIni": "0.1",
+        "stopCritIniSmall": "1.001",
+        "stopCritType": "kinEnergy",
+        "uFlowingThreshold": "0.1",
+        "gravAcc": "9.81",
+        "velMagMin": "1.e-6",
+        "rho": "100.",
+        "interpOption": "2",
+        "explicitFriction": "0",
+        "centeredPosition": "1",
+        "reprojMethodPosition": "2",
+        "reprojectionIterations": "5",
+        "thresholdProjection": "0.001",
+        "dissDam": "1",
+        "snowSlide": "1",
+        "wetSnow": "1",
+        "adaptSfcStopped": "0",
+        "adaptSfcDetrainment": "0",
+        "adaptSfcEntrainment": "0",
+    }
+
+    particles = {
+        "dt": 1.0,
+        "m": np.asarray([10.0, 10.0, 10.0], dtype=np.float64),
+        "idFixed": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "trajectoryLengthXY": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "trajectoryLengthXYCor": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "trajectoryLengthXYZ": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "x": np.asarray([0.0, 1.0, 2.0], dtype=np.float64),
+        "y": np.asarray([2.0, 3.0, 4.0], dtype=np.float64),
+        "z": np.asarray([1.0, 1.0, 1.0], dtype=np.float64),
+        "ux": np.asarray([1.0, 1.0, 1.0], dtype=np.float64),
+        "uy": np.asarray([1.0, 1.0, 1.0], dtype=np.float64),
+        "uz": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "uAcc": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "kineticEne": 0.0,
+        "peakKinEne": 0.0,
+        "peakForceSPH": 0.0,
+        "forceSPHIni": 0.0,
+        "nPart": 3,
+        "peakMassFlowing": 0.0,
+        "iterate": True,
+        "totalEnthalpy": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "velocityMag": np.asarray([1.0, 1.0, 1.0], dtype=np.float64),
+        "h": np.asarray([1.0, 1.0, 1.0], dtype=np.float64),
+        "ID": np.array([0, 1, 2], dtype=np.int64),
+        "stoppedParticles": {},
+        "indXDEM": np.array([0], dtype=np.int32),
+        "indYDEM": np.array([0], dtype=np.int32),
+    }
+    particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"])
+
+    demHeader = {}
+    demHeader["xllcenter"] = 0.0
+    demHeader["yllcenter"] = 0.0
+    demHeader["cellsize"] = 5.0
+    demHeader["nodata_value"] = -9999
+    demHeader["nrows"] = 10
+    demHeader["ncols"] = 10
+    dem = {"header": demHeader}
+    dem["rasterData"] = np.ones((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64)
+    dem["outOfDEM"] = np.zeros((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64).flatten()
+    dem["Nx"] = np.zeros((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64)
+    dem["Ny"] = np.zeros((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64)
+    dem["Nz"] = np.ones((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64)
+
+    force = {
+        "forceZ": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "forceFrict": np.asarray([10.0, 10.0, 10.0], dtype=np.float64),
+        "dM": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "dMDet": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+        "forceX": np.asarray([50.0, 50.0, 50.0], dtype=np.float64),
+        "forceY": np.asarray([50.0, 50.0, 50.0], dtype=np.float64),
+        "forceSPHX": np.asarray([50.0, 50.0, 50.0], dtype=np.float64),
+        "forceSPHY": np.asarray([50.0, 50.0, 50.0], dtype=np.float64),
+        "forceSPHZ": np.asarray([0.0, 0.0, 0.0], dtype=np.float64),
+    }
+    fields = {"FT": np.zeros((2, 2))}
+    wallLineDict = {
+        "dam": 0,
+        "nIterDam": 1,
+        "cellsCrossed": np.zeros((dem["header"]["ncols"] * dem["header"]["nrows"]), dtype=np.int64),
+    }
+    for key in ["x", "y", "z", "xCrown", "yCrown", "zCrown",
+                "xTangent", "yTangent", "zTangent"]:
+        wallLineDict[key] = np.ones((1)) * 1.0
+    for key in ["nPoints", "height", "slope", "restitutionCoefficient"]:
+        wallLineDict[key] = 0
+    dem["damLine"] = wallLineDict
+
+    particles = DFAfunC.updatePositionC(cfg["GENERAL"], particles, dem, force, fields, typeStop=0)
+
+    # Verify no circular reference was created
+    assert particles["stoppedParticles"] is not particles
+    assert "stoppedParticles" not in particles["stoppedParticles"]
+
+
 def test_computeTrajectoryAngle():
     # first compute travel angle for each particle
     # get parent Id in order to  get the first z position
@@ -619,6 +725,12 @@ def test_initializeBondsC():
     bondDist = particles["bondDist"]
     bondPart = particles["bondPart"]
     assert np.array_equal(bondStart, np.asarray([0, 2, 4, 6]))
+    assert particles["bondStart"].flags.owndata
+    assert particles["bondStart"].base is None
+    assert particles["bondPart"].flags.owndata
+    assert particles["bondPart"].base is None
+    assert particles["bondDist"].flags.owndata
+    assert particles["bondDist"].base is None
     for k in range(nPart):
         # loop on all bonded particles
         neighbors = list()
@@ -663,9 +775,18 @@ def test_removeBondsC():
     bondStart = particles["bondStart"]
     bondDist = particles["bondDist"]
     bondPart = particles["bondPart"]
-    assert np.array_equal(bondStart, np.asarray([0, 1, 1, 2]))
-    for k in range(nPart):
-        # loop on all bonded particles
+    # After removing particle 1 (old index), indices are remapped:
+    # old 0 -> new 0, removed, old 2 -> new 1
+    assert np.array_equal(bondStart, np.asarray([0, 1, 2]))
+    assert particles["bondStart"].flags.owndata
+    assert particles["bondStart"].base is None
+    assert particles["bondPart"].flags.owndata
+    assert particles["bondPart"].base is None
+    assert particles["bondDist"].flags.owndata
+    assert particles["bondDist"].base is None
+    # Only 2 particles remain; bondStart is new-indexed
+    nPartNew = nPart - nRemove
+    for k in range(nPartNew):
         neighbors = list()
         for ib in range(bondStart[k], bondStart[k + 1]):
             l = bondPart[ib]
@@ -673,15 +794,52 @@ def test_removeBondsC():
 
         neighbors.sort()
         if k == 0:
-            assert neighbors == [2]
+            assert neighbors == [1]
         if k == 1:
-            assert neighbors == []
-        if k == 2:
             assert neighbors == [0]
     bondDist.sort()
     assert np.array_equal(bondDist, np.asarray([1, 1]))
 
 
+def test_updateFieldsC_returns_owned_arrays():
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {"rho": "200.", "rhoEnt": "100.", "interpOption": "2"}
+    header = {"nrows": 5, "ncols": 5, "cellsize": 1}
+    dem = {"header": header, "areaRaster": 25 * np.ones((header["nrows"], header["ncols"]))}
+
+    particles = {
+        "m": np.array([100.0, 100.0, 100.0, 100.0]),
+        "dmDet": np.array([1.0, 1.0, 1.0, 1.0]),
+        "dmEnt": np.array([0.0, 0.0, 0.0, 0.0]),
+        "x": np.array([2.0, 1.0, 2.0, 1.0]),
+        "y": np.array([2.0, 2.0, 1.0, 1.0]),
+        "ux": np.array([10.0, 10.0, 10.0, 10.0]),
+        "uy": np.array([10.0, 10.0, 10.0, 10.0]),
+        "uz": np.array([10.0, 10.0, 10.0, 10.0]),
+        "trajectoryAngle": np.array([10.0, 10.0, 10.0, 10.0]),
+        "stoppedParticles": {"m": np.empty(0), "x": np.empty(0), "y": np.empty(0)},
+        "nPart": 4,
+    }
+
+    fields = {
+        "computeTA": False,
+        "computeKE": False,
+        "computeP": False,
+        "pfv": np.zeros((header["nrows"], header["ncols"])),
+        "ppr": np.zeros((header["nrows"], header["ncols"])),
+        "pft": np.zeros((header["nrows"], header["ncols"])),
+        "pta": np.zeros((header["nrows"], header["ncols"])),
+        "pke": np.zeros((header["nrows"], header["ncols"])),
+        "dmDet": np.zeros((header["nrows"], header["ncols"])),
+    }
+
+    particles, fields = DFAfunC.updateFieldsC(cfg["GENERAL"], particles, dem, fields)
+
+    for key in ["FM", "FV", "Vx", "Vy", "Vz", "FT", "pfv", "pft", "dmDet", "FTStop", "FTDet", "FTEnt"]:
+        assert fields[key].flags.owndata
+        assert fields[key].base is None
+
+
 def test_computeCohesionForceC():
     cfg = configparser.ConfigParser()
     cfg["GENERAL"] = {