diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx
index 568eb1268..8738b850e 100644
--- a/avaframe/com1DFA/DFAfunctionsCython.pyx
+++ b/avaframe/com1DFA/DFAfunctionsCython.pyx
@@ -1393,6 +1393,8 @@ def updateFieldsC(cfg, particles, dem, fields):
             PKE[j, i] = kineticEnergy[j, i]
 
       # thickness change due to detrainment, stopping and entrainment
+      # MassDetBilinear, MassStopBilinear <= 0, MassEntBilinear >= 0,
+      # FTDetBilinear, FTStopBilinear >= 0, FTEntBilinear <= 0
       FTDetBilinear[j, i] = - MassDetBilinear[j, i] / (areaRaster[j, i] * rho)  # / m * FTBilinear[j, i]
       FTStopBilinear[j, i] = - MassStopBilinear[j, i] / (areaRaster[j, i] * rho)  # / m * FTBilinear[j, i] 
       FTEntBilinear[j, i] = - MassEntBilinear[j, i] / (areaRaster[j, i] * rhoEnt)
@@ -1409,6 +1411,7 @@ def updateFieldsC(cfg, particles, dem, fields):
   fields['FTStop'] = np.asarray(FTStopBilinear).copy()
   fields['FTDet'] = np.asarray(FTDetBilinear).copy()
   fields['FTEnt'] = np.asarray(FTEntBilinear).copy()
+  fields['mStop'] = np.asarray(MassStopBilinear).copy()
 
   if computeP:
     fields['ppr'] = np.asarray(PP).copy()
diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py
index ec0cbb4c7..b4b8cceca 100644
--- a/avaframe/com1DFA/com1DFA.py
+++ b/avaframe/com1DFA/com1DFA.py
@@ -38,6 +38,7 @@
 import avaframe.com1DFA.DFAToolsCython as DFAtllsC
 import avaframe.com1DFA.damCom1DFA as damCom1DFA
 import avaframe.in2Trans.rasterUtils as IOf
+import avaframe.in2Trans.transformFields as transField
 import avaframe.in3Utils.fileHandlerUtils as fU
 from avaframe.in3Utils import cfgUtils
 import avaframe.out3Plot.outDebugPlots as debPlot
@@ -552,7 +553,7 @@ def setThickness(cfg, lineTh, typeTh):
     if cfg["GENERAL"].getboolean(thFlag):
         if cfg["INPUT"]["thFromIni"] != "" and typeTh in cfg["INPUT"]["thFromIni"]:
             lineTh["thicknessSource"] = ["ini file"] * len(lineTh["thickness"])
-        elif cfg["GENERAL"].getboolean("timeDependentRelease"):
+        elif cfg["GENERAL"].getboolean("timeDependentRelease") and typeTh == "relTh":
             lineTh["thicknessSource"] = ["csv file"] * len(lineTh["thickness"])
         else:
             lineTh["thicknessSource"] = ["shp file"] * len(lineTh["thickness"])
@@ -561,7 +562,7 @@ def setThickness(cfg, lineTh, typeTh):
 
     # set thickness value info read from ini file that has been updated from shp or ini previously
     for count, id in enumerate(lineTh["id"]):
-        if cfg["GENERAL"].getboolean("timeDependentRelease"):
+        if cfg["GENERAL"].getboolean("timeDependentRelease") and typeTh == "relTh":
             lineTh["thickness"][count] = float(lineTh["thickness"][count].item())
 
         elif cfg["GENERAL"].getboolean(thFlag):
@@ -1314,7 +1315,7 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     # initialize entrainment and resistance
     # get info of simType and whether or not to initialize resistance and entrainment
     simTypeActual = cfgGen["simTypeActual"]
-    entrMassRaster, entrEnthRaster, reportAreaInfo = initializeMassEnt(
+    entrMassRaster, entrEnthRaster, entrDepthRaster, reportAreaInfo = initializeMassEnt(
         dem,
         simTypeActual,
         inputSimLines["entLine"],
@@ -1326,6 +1327,12 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     # check if entrainment and release overlap
     entrMassRaster = geoTrans.checkOverlap(entrMassRaster, relRaster, "Entrainment", "Release", crop=True)
     entrEnthRaster = geoTrans.checkOverlap(entrEnthRaster, relRaster, "Entrainment", "Release", crop=True)
+
+    if cfg["GENERAL"].getboolean("entrainableDeposition"):
+        demNotErodableRaster = dem["rasterData"] - entrDepthRaster
+    else:
+        demNotErodableRaster = dem["rasterData"]
+
     # check for overlap with the secondary release area
     if secondaryReleaseInfo["flagSecondaryRelease"] == "Yes":
         for secIndex, secRelRaster in enumerate(secondaryReleaseInfo["rasterData"]):
@@ -1385,6 +1392,8 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     # surfacic entrainment mass available (unit kg/m²)
     fields["entrMassRaster"] = entrMassRaster
     fields["entrEnthRaster"] = entrEnthRaster
+    fields["entrDepth"] = entrDepthRaster
+    fields["demNotErodableRaster"] = demNotErodableRaster
     entreainableMass = np.nansum(fields["entrMassRaster"] * dem["areaRaster"])
     log.info("Mass available for entrainment: %.2f kg" % (entreainableMass))
 
@@ -1963,9 +1972,14 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin
     -------
     entrMassRaster : 2D numpy array
         raster of available mass for entrainment
+    entrEnthRaster: 2D numpy array
+        raster of available entrainment enthalpie
+    entrDepthDict["rasterData"]: 2D numpy array
+        raster of available depth for entrainment
     reportAreaInfo: dict
         simulation area information dictionary completed with entrainment area info
     """
+
     # read dem header
     header = dem["originalHeader"]
     ncols = header["ncols"]
@@ -1976,26 +1990,29 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin
         log.info("Entrainment area features: %s" % (entLine["Name"]))
         if entLine["initializedFrom"] == "shapefile":
             entLine = geoTrans.prepareArea(entLine, dem, thresholdPointInPoly, thList=entLine["thickness"])
-        entrMassRaster = entLine["rasterData"]
+        entrThicknessRaster = entLine["rasterData"]
         # ToDo: not used in samos but implemented
         # tempRaster = cfg['GENERAL'].getfloat('entTempRef') + (dem['rasterData'] - cfg['GENERAL'].getfloat('entMinZ'))
         # * cfg['GENERAL'].getfloat('entTempGrad')
         # entrEnthRaster = np.where(tempRaster < 0, tempRaster*cfg['GENERAL'].getfloat('cpIce'),
         #                           tempRaster*cfg['GENERAL'].getfloat('cpWtr')/cfg['GENERAL'].getfloat('hFusion'))
         entrEnthRaster = np.where(
-            entrMassRaster > 0,
+            entrThicknessRaster > 0,
             cfg["GENERAL"].getfloat("entTempRef") * cfg["GENERAL"].getfloat("cpIce"),
             0,
         )
         reportAreaInfo["entrainment"] = "Yes"
     else:
-        entrMassRaster = np.zeros((nrows, ncols))
+        entrThicknessRaster = np.zeros((nrows, ncols))
         entrEnthRaster = np.zeros((nrows, ncols))
         reportAreaInfo["entrainment"] = "No"
 
-    entrMassRaster = entrMassRaster * cfg["GENERAL"].getfloat("rhoEnt")
+    entrMassRaster = entrThicknessRaster * cfg["GENERAL"].getfloat("rhoEnt")
+    # convert entrainable thickness into entrainable depth
+    entrThDict = {"rasterData": entrThicknessRaster, "header": dem["header"]}
+    entrDepthDict, _, _ = transField.convertDepthThickness(entrThDict, dem, typeOfInput="thickness")
 
-    return entrMassRaster, entrEnthRaster, reportAreaInfo
+    return entrMassRaster, entrEnthRaster, entrDepthDict["rasterData"], reportAreaInfo
 
 
 def initializeResistance(cfg, dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly):
@@ -3499,6 +3516,20 @@ def prepareVarSimDict(standardCfg, inputSimFiles, variationDict, simNameExisting
                 cfgSim["INPUT"]["entThFile"] = pathToEnt
                 inputSimFiles["entResInfo"]["entRemeshed"] = remeshedEnt
             cfgSim["INPUT"]["entrainmentScenario"] = str(pathlib.Path("ENT", inputSimFiles["entFile"].name))
+            if cfgSim["GENERAL"]["entrainableDeposition"] == "True" and (
+                cfgSim["GENERAL"]["adaptSfcEntrainment"] == "0"
+                or (
+                    cfgSim["GENERAL"]["adaptSfcStopped"] == "0"
+                    and cfgSim["GENERAL"]["adaptSfcDetrainment"] == "0"
+                )
+            ):
+                # set entrainable Deposition to False if entrainment or adapt Sfc are switched off
+                message = f"When entrainable deposition is True, adaptSfcEntrainment and adaptSfcDeposition (or adaptSfcDetrainment) need to be 1."
+                log.error(message)
+                raise ValueError(message)
+        else:
+            cfgSim["GENERAL"]["entrainableDeposition"] = "False"
+            log.info("Deposition is not entrainable when it's not an ent sim type.")
 
         # add info about resistance file path to the cfg
         if "res" in row._asdict()["simTypeList"] and inputSimFiles["resFile"] is not None:
@@ -3907,7 +3938,15 @@ def adaptDEM(dem, fields, cfg):
     dem: dict
         dictionary with info on DEM data
     fields : dict
-        fields dictionary
+        fields dictionary containing:
+            FTDet (flow thickness change due to detrainment >= 0)
+            FTStop (flow thickness change due to stopping >= 0)
+            FTEnt (flow thickness change due to entrainment <= 0)
+            mStop (stopped mass <=0)
+            sfcChangeTotal (so far changed surface)
+            entrDepth (entrainable depth that is on top of the not erodable topography)
+            entrMassRaster (entrainable mass)
+            demNotErodableRaster (topography that is not erodable)
     cfg: dict
         configuration settings
 
@@ -3916,17 +3955,23 @@ def adaptDEM(dem, fields, cfg):
     dem: dict
         dictionary with info on DEM data containing adapted topography
     fields : dict
-        fields dictionary containing adapted DEM
+        fields dictionary containing adapted DEM, adapted entrainable mass and depth raster
     """
 
     ZDEM = dem["rasterData"].copy()
-    FTDet = fields["FTDet"]
-    FTStop = fields["FTStop"]
-    FTEnt = fields["FTEnt"]
+    FTDet = fields["FTDet"]  # >= 0
+    FTStop = fields["FTStop"]  # >= 0
+    FTEnt = fields["FTEnt"]  # <= 0
+    mStop = fields["mStop"]  # <= 0
+
     sfcChangeTotal = fields["sfcChangeTotal"]
     sfcChange = np.zeros_like(FTDet)
     ZDEMadapt = ZDEM
 
+    depthStop = np.zeros_like(ZDEM)
+    depthDet = np.zeros_like(ZDEM)
+    depthEnt = np.zeros_like(ZDEM)
+
     _, _, NzNormed = DFAtls.normalize(dem["Nx"].copy(), dem["Ny"].copy(), dem["Nz"].copy())
 
     if cfg.getboolean("adaptSfcStopped"):
@@ -3944,6 +3989,14 @@ def adaptDEM(dem, fields, cfg):
         depthEnt = FTEnt / NzNormed
         ZDEMadapt += depthEnt
         sfcChange += depthEnt
+    if cfg.getboolean("entrainableDeposition"):
+        hErodable = fields["entrDepth"]
+        hErodable += depthStop
+        hErodable += depthDet
+        hErodable += depthEnt
+        ZDEMadapt = fields["demNotErodableRaster"] + hErodable
+        fields["entrDepth"] = hErodable
+        fields["entrMassRaster"] -= mStop
 
     dem["rasterData"] = ZDEMadapt
     fields["demAdapted"] = ZDEMadapt
diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini
index 5669d0428..884f826a9 100644
--- a/avaframe/com1DFA/com1DFACfg.ini
+++ b/avaframe/com1DFA/com1DFACfg.ini
@@ -11,7 +11,7 @@ simTypeList = available
 modelType = dfa
 
 #+++++++++++++ Output++++++++++++
-# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, dmDet, sfcChange, demAdapted, timeInfo, particles) - separated by |
+# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, dmDet, sfcChange, demAdapted, timeInfo, entrDepth, particles) - separated by |
 resType = ppr|pft|pfv|timeInfo
 # saving time step, i.e.  time in seconds (last time step is always saved; initial time step only if explicitly specified)
 # option 1: give an interval with start:interval in seconds (tStep = 0:5 - this will save desired results every 5 seconds for the full simulation)
@@ -489,6 +489,9 @@ entDefResistance = 0
 adaptSfcStopped = 0
 adaptSfcDetrainment = 0
 adaptSfcEntrainment = 0
+# deposition is only entrainable if it is an entrainment simulation and if adaptSfcStopped = 1 and adaptSfcEntrainment = 1
+# then mass that is stopped can be entrained (also if the entrainment area is not in this stopped area)
+entrainableDeposition = False
 # only adapt topography if changing height in at least one cell is > thresholdAdaptSfc [m]
 thresholdAdaptSfc = 0.1
 # use the adapted topography as background in the report plots
diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py
index 6b536284d..2e0b2e512 100644
--- a/avaframe/com1DFA/deriveParameterSet.py
+++ b/avaframe/com1DFA/deriveParameterSet.py
@@ -154,6 +154,7 @@ def checkResType(fullCfg, section, key, value):
             "sfcChange",
             "demAdapted",
             "timeInfo",
+            "entrDepth",
         ]
         message = "The parameter % s is not a valid resType. It will not be saved"
         newResType = []
@@ -285,7 +286,7 @@ def getThicknessValue(cfg, inputSimFiles, fName, thType):
 
     # if thickness should be read from shape file
     if cfg["GENERAL"].getboolean(thFlag):
-        if cfg["GENERAL"].getboolean("timeDependentRelease"):
+        if cfg["GENERAL"].getboolean("timeDependentRelease") and thType == "relTh":
             cfg["INPUT"]["relThInfo"] = "timeDependent"
         # if at least one but not all features in a shapefile have a thickness value - error
         elif ("None" in thicknessList) and thType != "entTh":
@@ -902,8 +903,11 @@ def appendThicknessToCfg(cfg):
             for count, id in enumerate(idList):
                 thNameId = thType + id
                 if thNameId in cfg["GENERAL"].keys():
-                    log.info("Thickness value for %s already set in initial config file, \
-                              read from there not from shp file" % thNameId)
+                    log.info(
+                        "Thickness value for %s already set in initial config file, \
+                              read from there not from shp file"
+                        % thNameId
+                    )
                 else:
                     cfgGen[thNameId] = str(float(thicknessList[count]))
 
diff --git a/avaframe/in2Trans/transformFields.py b/avaframe/in2Trans/transformFields.py
index de711aa6f..2f15775de 100644
--- a/avaframe/in2Trans/transformFields.py
+++ b/avaframe/in2Trans/transformFields.py
@@ -1,5 +1,5 @@
 """
-    ASCII file reader and handler
+ASCII file reader and handler
 
 """
 
@@ -16,14 +16,14 @@
 log = logging.getLogger(__name__)
 
 
-def convertDepthThickness(inputDict, demDict, typeOfInput='depth'):
+def convertDepthThickness(inputDict, demDictOrig, typeOfInput="depth"):
     """convert depthField to thicknessField or the other way around depending on type, using a DEM to compute the slope
     angle required for transformation
     also writes field to a new directory called transformed where depthFile/thicknessFile is located
 
     Parameters
     -----------
-    demDict: dict
+    demDictOrig: dict
         dictionary with dem header and rasterData (numpy nd array of z values)
     inputDict: dict
         dictionary with depthField/thicknessField header and rasterData (numpy nd array of depth values)
@@ -37,6 +37,7 @@ def convertDepthThickness(inputDict, demDict, typeOfInput='depth'):
 
     """
     # get normal vector of the grid mesh
+    demDict = {"header": demDictOrig["header"].copy(), "rasterData": demDictOrig["rasterData"]}
     demDict = gT.getNormalMesh(demDict)
     _, _, NzNormed = DFAtls.normalize(demDict["Nx"], demDict["Ny"], demDict["Nz"])
 
@@ -45,41 +46,41 @@ def convertDepthThickness(inputDict, demDict, typeOfInput='depth'):
     demDict["header"]["nodata_value"] = inputDict["header"]["nodata_value"]
 
     # multiply depth with cos(slopeAngle)
-    if typeOfInput == 'depth':
+    if typeOfInput == "depth":
         outField = inputRasterNew * NzNormed
-    elif typeOfInput == 'thickness':
+    elif typeOfInput == "thickness":
         outField = inputRasterNew / NzNormed
     else:
-        message = 'Type for thickness/depth conversion is: %s - not valid' % typeOfInput
+        message = "Type for thickness/depth conversion is: %s - not valid" % typeOfInput
         log.error(message)
         raise AssertionError(message)
 
     slopeAngleField = np.rad2deg(np.arccos(NzNormed))
 
     # create thickness dict
-    outputDict = {'header': demDict['header'], 'rasterData': outField}
+    outputDict = {"header": demDict["header"], "rasterData": outField}
 
     return outputDict, inputRasterNew, slopeAngleField
 
 
-def fetchPointValuesFromField(dataDF, xyPoints, resType, interpMethod='bilinear'):
-    """ derive field values at xyPoints using a interpMethod (options: nearest and bilinear)
+def fetchPointValuesFromField(dataDF, xyPoints, resType, interpMethod="bilinear"):
+    """derive field values at xyPoints using a interpMethod (options: nearest and bilinear)
 
-        Parameters
-        -----------
-        dataDF: pandas dataFrame
-            dataframe with info on simulations (including result file paths,configuration)
-        xyPoints: dict
-            dictionary with keys x, y and point of interest coordinates
-        resType: str
-            name of result type (pfv, pft, ...)
-        interpMethod: str
-            interpolation method to derive values of field at xyPoints
+    Parameters
+    -----------
+    dataDF: pandas dataFrame
+        dataframe with info on simulations (including result file paths,configuration)
+    xyPoints: dict
+        dictionary with keys x, y and point of interest coordinates
+    resType: str
+        name of result type (pfv, pft, ...)
+    interpMethod: str
+        interpolation method to derive values of field at xyPoints
 
-        Returns
-        ---------
-        dataDF: pandas dataFrame
-            updated pandas dataFrame with new column pointValues_resType providing value of resType at xyPoint
+    Returns
+    ---------
+    dataDF: pandas dataFrame
+        updated pandas dataFrame with new column pointValues_resType providing value of resType at xyPoint
 
     """
 
@@ -88,9 +89,9 @@ def fetchPointValuesFromField(dataDF, xyPoints, resType, interpMethod='bilinear'
         # read field
         field = IOf.readRaster(row[resType])
 
-        value, _ = gT.projectOnRaster(field, xyPoints, interp=interpMethod, inData="rasterData", outData="value")
-        dataDF.loc[index, ('pointValues_%s' % resType)] = value['value'][0]
+        value, _ = gT.projectOnRaster(
+            field, xyPoints, interp=interpMethod, inData="rasterData", outData="value"
+        )
+        dataDF.loc[index, ("pointValues_%s" % resType)] = value["value"][0]
 
     return dataDF
-
-
diff --git a/avaframe/in3Utils/geoTrans.py b/avaframe/in3Utils/geoTrans.py
index fb3ca157a..12e6984b8 100644
--- a/avaframe/in3Utils/geoTrans.py
+++ b/avaframe/in3Utils/geoTrans.py
@@ -1,5 +1,4 @@
-""" Opperations and transformations of rasters and lines
-"""
+"""Opperations and transformations of rasters and lines"""
 
 import logging
 import math
@@ -445,7 +444,7 @@ def remeshRaster(rasterFile, cfgSim, typeIndicator="DEM", onlySearch=False, lega
     fU.makeADir(pathToRaster)
 
     outFile = pathToRaster / (
-            "%s_remeshed%s%.2f" % (rasterFile.stem, typeIndicator, remeshedRaster["header"]["cellsize"])
+        "%s_remeshed%s%.2f" % (rasterFile.stem, typeIndicator, remeshedRaster["header"]["cellsize"])
     )
 
     writtenFile = rU.writeResultToRaster(
@@ -497,8 +496,8 @@ def searchRemeshedRaster(rasterName, cfgSim, typeIndicator="DEM"):
         for rasterF in rasterFiles:
             headerRaster = rU.readRasterHeader(rasterF)
             if (
-                    abs(meshCellSize - headerRaster["cellsize"]) < meshCellSizeThreshold
-                    and rasterName in rasterF.stem
+                abs(meshCellSize - headerRaster["cellsize"]) < meshCellSizeThreshold
+                and rasterName in rasterF.stem
             ):
                 log.info(
                     "Remeshed Raster found: %s cellSize: %.5f" % (rasterF.name, headerRaster["cellsize"])
@@ -535,8 +534,8 @@ def computeS(avaPath):
     ycoord = avaPath["y"]
     n = np.size(xcoord)
     # compute s
-    dxs = xcoord[1:n] - xcoord[0: n - 1]
-    dys = ycoord[1:n] - ycoord[0: n - 1]
+    dxs = xcoord[1:n] - xcoord[0 : n - 1]
+    dys = ycoord[1:n] - ycoord[0 : n - 1]
     # deduce the distance in s direction
     ds2 = dxs * dxs + dys * dys
     ds = np.sqrt(ds2)
@@ -1786,30 +1785,30 @@ def getNormalMesh(dem, num=4):
         # filling the inside of the matrix
         # normal calculation with 4 triangles
         # (Zl - Zr) / csz
-        Nx[1: n - 1, 1: m - 1] = (z[1: n - 1, 0: m - 2] - z[1: n - 1, 2:m]) / csz
+        Nx[1 : n - 1, 1 : m - 1] = (z[1 : n - 1, 0 : m - 2] - z[1 : n - 1, 2:m]) / csz
         # (Zd - Zu) * csz
-        Ny[1: n - 1, 1: m - 1] = (z[0: n - 2, 1: m - 1] - z[2:n, 1: m - 1]) / csz
+        Ny[1 : n - 1, 1 : m - 1] = (z[0 : n - 2, 1 : m - 1] - z[2:n, 1 : m - 1]) / csz
         Nz = 2 * Nz
         # filling the first col of the matrix
         # -2*(Zr - Zp) / csz
-        Nx[1: n - 1, 0] = -2 * (z[1: n - 1, 1] - z[1: n - 1, 0]) / csz
+        Nx[1 : n - 1, 0] = -2 * (z[1 : n - 1, 1] - z[1 : n - 1, 0]) / csz
         # (Zd - Zu) / csz
-        Ny[1: n - 1, 0] = (z[0: n - 2, 0] - z[2:n, 0]) / csz
+        Ny[1 : n - 1, 0] = (z[0 : n - 2, 0] - z[2:n, 0]) / csz
         # filling the last col of the matrix
         # 2*(Zl - Zp) / csz
-        Nx[1: n - 1, m - 1] = 2 * (z[1: n - 1, m - 2] - z[1: n - 1, m - 1]) / csz
+        Nx[1 : n - 1, m - 1] = 2 * (z[1 : n - 1, m - 2] - z[1 : n - 1, m - 1]) / csz
         # (Zd - Zu) / csz
-        Ny[1: n - 1, m - 1] = (z[0: n - 2, m - 1] - z[2:n, m - 1]) / csz
+        Ny[1 : n - 1, m - 1] = (z[0 : n - 2, m - 1] - z[2:n, m - 1]) / csz
         # filling the first row of the matrix
         # (Zl - Zr) / csz
-        Nx[0, 1: m - 1] = (z[0, 0: m - 2] - z[0, 2:m]) / csz
+        Nx[0, 1 : m - 1] = (z[0, 0 : m - 2] - z[0, 2:m]) / csz
         # -2*(Zu - Zp) / csz
-        Ny[0, 1: m - 1] = -2 * (z[1, 1: m - 1] - z[0, 1: m - 1]) / csz
+        Ny[0, 1 : m - 1] = -2 * (z[1, 1 : m - 1] - z[0, 1 : m - 1]) / csz
         # filling the last row of the matrix
         # (Zl - Zr) / csz
-        Nx[n - 1, 1: m - 1] = (z[n - 1, 0: m - 2] - z[n - 1, 2:m]) / csz
+        Nx[n - 1, 1 : m - 1] = (z[n - 1, 0 : m - 2] - z[n - 1, 2:m]) / csz
         # 2*(Zd - Zp) / csz
-        Ny[n - 1, 1: m - 1] = 2 * (z[n - 2, 1: m - 1] - z[n - 1, 1: m - 1]) / csz
+        Ny[n - 1, 1 : m - 1] = 2 * (z[n - 2, 1 : m - 1] - z[n - 1, 1 : m - 1]) / csz
         # filling the corners of the matrix
         Nx[0, 0] = -(z[0, 1] - z[0, 0]) / csz
         Ny[0, 0] = -(z[1, 0] - z[0, 0]) / csz
@@ -1828,56 +1827,54 @@ def getNormalMesh(dem, num=4):
         # filling the inside of the matrix
         # normal calculation with 6 triangles
         # (2*(Zl - Zr) - Zur + Zdl + Zu - Zd) / csz
-        Nx[1: n - 1, 1: m - 1] = (
-                                         2 * (z[1: n - 1, 0: m - 2] - z[1: n - 1, 2:m])
-                                         - z[2:n, 2:m]
-                                         + z[0: n - 2, 0: m - 2]
-                                         + z[2:n, 1: m - 1]
-                                         - z[0: n - 2, 1: m - 1]
-                                 ) / csz
+        Nx[1 : n - 1, 1 : m - 1] = (
+            2 * (z[1 : n - 1, 0 : m - 2] - z[1 : n - 1, 2:m])
+            - z[2:n, 2:m]
+            + z[0 : n - 2, 0 : m - 2]
+            + z[2:n, 1 : m - 1]
+            - z[0 : n - 2, 1 : m - 1]
+        ) / csz
         # (2*(Zd - Zu) - Zur + Zdl - Zl + Zr) / csz
-        Ny[1: n - 1, 1: m - 1] = (
-                                         2 * (z[0: n - 2, 1: m - 1] - z[2:n, 1: m - 1])
-                                         - z[2:n, 2:m]
-                                         + z[0: n - 2, 0: m - 2]
-                                         - z[1: n - 1, 0: m - 2]
-                                         + z[1: n - 1, 2:m]
-                                 ) / csz
+        Ny[1 : n - 1, 1 : m - 1] = (
+            2 * (z[0 : n - 2, 1 : m - 1] - z[2:n, 1 : m - 1])
+            - z[2:n, 2:m]
+            + z[0 : n - 2, 0 : m - 2]
+            - z[1 : n - 1, 0 : m - 2]
+            + z[1 : n - 1, 2:m]
+        ) / csz
         Nz = 6 * Nz
         # filling the first col of the matrix
         # (- 2*(Zr - Zp) + Zu - Zur ) / csz
-        Nx[1: n - 1, 0] = (-2 * (z[1: n - 1, 1] - z[1: n - 1, 0]) + z[2:n, 0] - z[2:n, 1]) / csz
+        Nx[1 : n - 1, 0] = (-2 * (z[1 : n - 1, 1] - z[1 : n - 1, 0]) + z[2:n, 0] - z[2:n, 1]) / csz
         # (Zd - Zu + Zr - Zur) / csz
-        Ny[1: n - 1, 0] = (z[0: n - 2, 0] - z[2:n, 0] + z[1: n - 1, 1] - z[2:n, 1]) / csz
-        Nz[1: n - 1, 0] = 3
+        Ny[1 : n - 1, 0] = (z[0 : n - 2, 0] - z[2:n, 0] + z[1 : n - 1, 1] - z[2:n, 1]) / csz
+        Nz[1 : n - 1, 0] = 3
         # filling the last col of the matrix
         # (2*(Zl - Zp) + Zdl - Zd) / csz
-        Nx[1: n - 1, m - 1] = (
-                                      2 * (z[1: n - 1, m - 2] - z[1: n - 1, m - 1]) + z[0: n - 2, m - 2] - z[0: n - 2,
-                                                                                                           m - 1]
-                              ) / csz
+        Nx[1 : n - 1, m - 1] = (
+            2 * (z[1 : n - 1, m - 2] - z[1 : n - 1, m - 1]) + z[0 : n - 2, m - 2] - z[0 : n - 2, m - 1]
+        ) / csz
         # (Zd - Zu + Zdl - Zl) / csz
-        Ny[1: n - 1, m - 1] = (
-                                      z[0: n - 2, m - 1] - z[2:n, m - 1] + z[0: n - 2, m - 2] - z[1: n - 1, m - 2]
-                              ) / csz
-        Nz[1: n - 1, m - 1] = 3
+        Ny[1 : n - 1, m - 1] = (
+            z[0 : n - 2, m - 1] - z[2:n, m - 1] + z[0 : n - 2, m - 2] - z[1 : n - 1, m - 2]
+        ) / csz
+        Nz[1 : n - 1, m - 1] = 3
         # filling the first row of the matrix
         # (Zl - Zr + Zu - Zur) / csz
-        Nx[0, 1: m - 1] = (z[0, 0: m - 2] - z[0, 2:m] + z[1, 1: m - 1] - z[1, 2:m]) / csz
+        Nx[0, 1 : m - 1] = (z[0, 0 : m - 2] - z[0, 2:m] + z[1, 1 : m - 1] - z[1, 2:m]) / csz
         # (-2*(Zu - Zp) + Zr - Zur) / csz
-        Ny[0, 1: m - 1] = (-2 * (z[1, 1: m - 1] - z[0, 1: m - 1]) + z[0, 2:m] - z[1, 2:m]) / csz
-        Nz[0, 1: m - 1] = 3
+        Ny[0, 1 : m - 1] = (-2 * (z[1, 1 : m - 1] - z[0, 1 : m - 1]) + z[0, 2:m] - z[1, 2:m]) / csz
+        Nz[0, 1 : m - 1] = 3
         # filling the last row of the matrix
         # (Zl - Zr + Zdl - Zd) / csz
-        Nx[n - 1, 1: m - 1] = (
-                                      z[n - 1, 0: m - 2] - z[n - 1, 2:m] + z[n - 2, 0: m - 2] - z[n - 2, 1: m - 1]
-                              ) / csz
+        Nx[n - 1, 1 : m - 1] = (
+            z[n - 1, 0 : m - 2] - z[n - 1, 2:m] + z[n - 2, 0 : m - 2] - z[n - 2, 1 : m - 1]
+        ) / csz
         # (2*(Zd - Zp) + Zdl - Zl) / csz
-        Ny[n - 1, 1: m - 1] = (
-                                      2 * (z[n - 2, 1: m - 1] - z[n - 1, 1: m - 1]) + z[n - 2, 0: m - 2] - z[n - 1,
-                                                                                                           0: m - 2]
-                              ) / csz
-        Nz[n - 1, 1: m - 1] = 3
+        Ny[n - 1, 1 : m - 1] = (
+            2 * (z[n - 2, 1 : m - 1] - z[n - 1, 1 : m - 1]) + z[n - 2, 0 : m - 2] - z[n - 1, 0 : m - 2]
+        ) / csz
+        Nz[n - 1, 1 : m - 1] = 3
         # filling the corners of the matrix
         Nx[0, 0] = (z[1, 0] - z[1, 1] - (z[0, 1] - z[0, 0])) / csz
         Ny[0, 0] = (z[0, 1] - z[1, 1] - (z[1, 0] - z[0, 0])) / csz
@@ -1896,74 +1893,74 @@ def getNormalMesh(dem, num=4):
         # filling the inside of the matrix
         # normal calculation with 8 triangles
         # (2*(Zl - Zr) + Zul - Zur + Zdl - Zdr) / csz
-        Nx[1: n - 1, 1: m - 1] = (
-                                         2 * (z[1: n - 1, 0: m - 2] - z[1: n - 1, 2:m])
-                                         + z[2:n, 0: m - 2]
-                                         - z[2:n, 2:m]
-                                         + z[0: n - 2, 0: m - 2]
-                                         - z[0: n - 2, 2:m]
-                                 ) / csz
+        Nx[1 : n - 1, 1 : m - 1] = (
+            2 * (z[1 : n - 1, 0 : m - 2] - z[1 : n - 1, 2:m])
+            + z[2:n, 0 : m - 2]
+            - z[2:n, 2:m]
+            + z[0 : n - 2, 0 : m - 2]
+            - z[0 : n - 2, 2:m]
+        ) / csz
         # (2*(Zd - Zu) - Zul - Zur + Zdl + Zdr) / csz
-        Ny[1: n - 1, 1: m - 1] = (
-                                         2 * (z[0: n - 2, 1: m - 1] - z[2:n, 1: m - 1])
-                                         - z[2:n, 0: m - 2]
-                                         - z[2:n, 2:m]
-                                         + z[0: n - 2, 0: m - 2]
-                                         + z[0: n - 2, 2:m]
-                                 ) / csz
+        Ny[1 : n - 1, 1 : m - 1] = (
+            2 * (z[0 : n - 2, 1 : m - 1] - z[2:n, 1 : m - 1])
+            - z[2:n, 0 : m - 2]
+            - z[2:n, 2:m]
+            + z[0 : n - 2, 0 : m - 2]
+            + z[0 : n - 2, 2:m]
+        ) / csz
         Nz = 8 * Nz
         # filling the first col of the matrix
         # (- 2*(Zr - Zp) + Zu - Zur + Zd - Zdr) / csz
-        Nx[1: n - 1, 0] = (
-                                  -2 * (z[1: n - 1, 1] - z[1: n - 1, 0])
-                                  + z[2:n, 0]
-                                  - z[2:n, 1]
-                                  + z[0: n - 2, 0]
-                                  - z[0: n - 2, 1]
-                          ) / csz
+        Nx[1 : n - 1, 0] = (
+            -2 * (z[1 : n - 1, 1] - z[1 : n - 1, 0])
+            + z[2:n, 0]
+            - z[2:n, 1]
+            + z[0 : n - 2, 0]
+            - z[0 : n - 2, 1]
+        ) / csz
         # (Zd - Zu + Zdr - Zur) / csz
-        Ny[1: n - 1, 0] = (z[0: n - 2, 0] - z[2:n, 0] + z[0: n - 2, 1] - z[2:n, 1]) / csz
-        Nz[1: n - 1, 0] = 4
+        Ny[1 : n - 1, 0] = (z[0 : n - 2, 0] - z[2:n, 0] + z[0 : n - 2, 1] - z[2:n, 1]) / csz
+        Nz[1 : n - 1, 0] = 4
         # filling the last col of the matrix
         # (2*(Zl - Zp) + Zdl - Zd + Zul - Zu) / csz
-        Nx[1: n - 1, m - 1] = (
-                                      2 * (z[1: n - 1, m - 2] - z[1: n - 1, m - 1])
-                                      + z[0: n - 2, m - 2]
-                                      - z[0: n - 2, m - 1]
-                                      + z[2:n, m - 2]
-                                      - z[2:n, m - 1]
-                              ) / csz
+        Nx[1 : n - 1, m - 1] = (
+            2 * (z[1 : n - 1, m - 2] - z[1 : n - 1, m - 1])
+            + z[0 : n - 2, m - 2]
+            - z[0 : n - 2, m - 1]
+            + z[2:n, m - 2]
+            - z[2:n, m - 1]
+        ) / csz
         # (Zd - Zu + Zdl - Zul) / csz
-        Ny[1: n - 1, m - 1] = (
-                                      z[0: n - 2, m - 1] - z[2:n, m - 1] + z[0: n - 2, m - 2] - z[2:n, m - 2]
-                              ) / csz
-        Nz[1: n - 1, m - 1] = 4
+        Ny[1 : n - 1, m - 1] = (
+            z[0 : n - 2, m - 1] - z[2:n, m - 1] + z[0 : n - 2, m - 2] - z[2:n, m - 2]
+        ) / csz
+        Nz[1 : n - 1, m - 1] = 4
         # filling the first row of the matrix
         # (Zl - Zr + Zul - Zur) / csz
-        Nx[0, 1: m - 1] = (z[0, 0: m - 2] - z[0, 2:m] + z[1, 0: m - 2] - z[1, 2:m]) / csz
+        Nx[0, 1 : m - 1] = (z[0, 0 : m - 2] - z[0, 2:m] + z[1, 0 : m - 2] - z[1, 2:m]) / csz
         # (-2*(Zu - Zp) + Zr - Zur + Zl - Zul) / csz
-        Ny[0, 1: m - 1] = (
-                                  -2 * (z[1, 1: m - 1] - z[0, 1: m - 1])
-                                  + z[0, 2:m]
-                                  - z[1, 2:m]
-                                  + z[0, 0: m - 2]
-                                  - z[1, 0: m - 2]
-                          ) / csz
-        Nz[0, 1: m - 1] = 4
+        Ny[0, 1 : m - 1] = (
+            -2 * (z[1, 1 : m - 1] - z[0, 1 : m - 1])
+            + z[0, 2:m]
+            - z[1, 2:m]
+            + z[0, 0 : m - 2]
+            - z[1, 0 : m - 2]
+        ) / csz
+        Nz[0, 1 : m - 1] = 4
         # filling the last row of the matrix
         # (Zl - Zr + Zdl - Zdr) / csz
-        Nx[n - 1, 1: m - 1] = (
-                                      z[n - 1, 0: m - 2] - z[n - 1, 2:m] + z[n - 2, 0: m - 2] - z[n - 2, 2:m]
-                              ) / csz
+        Nx[n - 1, 1 : m - 1] = (
+            z[n - 1, 0 : m - 2] - z[n - 1, 2:m] + z[n - 2, 0 : m - 2] - z[n - 2, 2:m]
+        ) / csz
         # (2*(Zd - Zp) + Zdl - Zl + Zdr - Zr) / csz
-        Ny[n - 1, 1: m - 1] = (
-                                      2 * (z[n - 2, 1: m - 1] - z[n - 1, 1: m - 1])
-                                      + z[n - 2, 0: m - 2]
-                                      - z[n - 1, 0: m - 2]
-                                      + z[n - 2, 2:m]
-                                      - z[n - 1, 2:m]
-                              ) / csz
-        Nz[n - 1, 1: m - 1] = 4
+        Ny[n - 1, 1 : m - 1] = (
+            2 * (z[n - 2, 1 : m - 1] - z[n - 1, 1 : m - 1])
+            + z[n - 2, 0 : m - 2]
+            - z[n - 1, 0 : m - 2]
+            + z[n - 2, 2:m]
+            - z[n - 1, 2:m]
+        ) / csz
+        Nz[n - 1, 1 : m - 1] = 4
         # filling the corners of the matrix
         Nx[0, 0] = (z[1, 0] - z[1, 1] - (z[0, 1] - z[0, 0])) / csz
         Ny[0, 0] = (z[0, 1] - z[1, 1] - (z[1, 0] - z[0, 0])) / csz
@@ -1986,19 +1983,19 @@ def getNormalMesh(dem, num=4):
         n2, m2 = np.shape(z2)
 
         Nx = (
-                -((z2[0: n2 - 1, 1:m2] - z2[1:n2, 0: m2 - 1]) + (z2[1:n2, 1:m2] - z2[0: n2 - 1, 0: m2 - 1]))
-                * csz
+            -((z2[0 : n2 - 1, 1:m2] - z2[1:n2, 0 : m2 - 1]) + (z2[1:n2, 1:m2] - z2[0 : n2 - 1, 0 : m2 - 1]))
+            * csz
         )
         Ny = (
-                -((z2[1:n2, 1:m2] - z2[0: n2 - 1, 0: m2 - 1]) - (z2[0: n2 - 1, 1:m2] - z2[1:n2, 0: m2 - 1]))
-                * csz
+            -((z2[1:n2, 1:m2] - z2[0 : n2 - 1, 0 : m2 - 1]) - (z2[0 : n2 - 1, 1:m2] - z2[1:n2, 0 : m2 - 1]))
+            * csz
         )
         Nz = 2 * Nz * csz * csz
 
         # Nx = - (z2[0:n2-1, 1:m2] - z2[0:n2-1, 0:m2-1]) / csz
         # Ny = - (z2[1:n2, 0:m2-1] - z2[0:n2-1, 0:m2-1]) / csz
-        Ny[n - 1, 0: m - 1] = -Ny[n - 1, 0: m - 1]
-        Nx[0: n - 1, m - 1] = -Nx[0: n - 1, m - 1]
+        Ny[n - 1, 0 : m - 1] = -Ny[n - 1, 0 : m - 1]
+        Nx[0 : n - 1, m - 1] = -Nx[0 : n - 1, m - 1]
         Ny[n - 1, m - 1] = -Ny[n - 1, m - 1]
         Nx[n - 1, m - 1] = -Nx[n - 1, m - 1]
         # TODO, Try to replicate samosAT notmal computation
@@ -2132,4 +2129,3 @@ def checkDBOverlap(DBXl, DBXr, DBYl, DBYr):
     if not DBrLine.is_simple or not DBlLine.is_simple:
         message = "Domain transformation for given path_aimec - curvature of provided line leads to folding"
         log.warning(message)
-
diff --git a/avaframe/out3Plot/plotUtils.py b/avaframe/out3Plot/plotUtils.py
index 27664f165..c239b6631 100644
--- a/avaframe/out3Plot/plotUtils.py
+++ b/avaframe/out3Plot/plotUtils.py
@@ -216,6 +216,18 @@
 ]
 cmapSfcChange = copy.copy(cmapCrameri.nuuk.reversed())
 
+levEntrDepth = list(fU.splitIniValueToArraySteps(cfgPlotUtils["entrDepthLevels"]))
+colorsEntrDepth = [
+    "#fefeb2",
+    "#d2d184",
+    "#bab98d",
+    "#a0a598",
+    "#6f878d",
+    "#386982",
+    "#05598c",
+]
+cmapEntrDepth = copy.copy(cmapCrameri.nuuk.reversed())
+
 # colormap for timeInfo
 levtimeInfo = list(fU.splitIniValueToArraySteps(cfgPlotUtils["timeInfoColorLevels"]))
 # lipari reversed color map
@@ -264,6 +276,12 @@
 
 cmapSfcChange = {"cmap": cmapSfcChange, "colors": colorsSfcC, "levels": levSfcC}
 
+cmapEntrDepth = {
+    "cmap": cmapEntrDepth,
+    "colors": colorsEntrDepth,
+    "levels": levEntrDepth,
+}
+
 cmapTime = {"cmap": cmaptimeInfo, "colors": colorstimInfo, "levels": levtimeInfo}
 
 # for zdelta
@@ -298,6 +316,7 @@
     "dmDet": cmapDmDet,
     "demAdapted": cmapGreys,
     "sfcChange": cmapSfcChange,
+    "entrDepth": cmapEntrDepth,
     "timeInfo": cmapTime,
 }
 
diff --git a/avaframe/out3Plot/plotUtilsCfg.ini b/avaframe/out3Plot/plotUtilsCfg.ini
index 277bbfdd5..a75ba3753 100644
--- a/avaframe/out3Plot/plotUtilsCfg.ini
+++ b/avaframe/out3Plot/plotUtilsCfg.ini
@@ -90,6 +90,7 @@ unitzdelta = m
 unitdmDet = kg
 unitsfcChange = m
 unitdemAdapted = m
+unitentrDepth = m
 unittimeInfo = s
 # threshold levels
 elevMaxppr = 100
@@ -117,6 +118,8 @@ travelAngleColorLevels = 28|29|30|31|32|33|34
 probaColorLevels = 0|0.25|0.50|0.75|1.
 # for thicknessChange
 surfaceChangeLevels = -1|-0.3|-0.1|0|0.1|0.3|1
+# for erodable depth
+entrDepthLevels = 0|0.5|1.0|1.5|2.0|2.5|3.0
 # for time Info
 timeInfoColorLevels = 10|20|30|40|50|60|70|80|90|100|110|120|130|140|150|160|170|180|200
 # contour levels (when adding contour lines on a plot)
diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py
index 17f1cd8cd..1661f8446 100644
--- a/avaframe/tests/test_com1DFA.py
+++ b/avaframe/tests/test_com1DFA.py
@@ -25,6 +25,7 @@
 import avaframe.com1DFA.DFAtools as DFAtls
 import avaframe.com1DFA.particleInitialisation as pI
 
+
 def test_prepareInputData(tmp_path):
     """test preparing input data"""
 
@@ -862,7 +863,7 @@ def test_initializeMassEnt():
     thresholdPointInPoly = 0.001
 
     # call function to be tested
-    entrMassRaster, entrEnthRaster, reportAreaInfo = com1DFA.initializeMassEnt(
+    entrMassRaster, entrEnthRaster, entrDepthRaster, reportAreaInfo = com1DFA.initializeMassEnt(
         dem,
         simTypeActual,
         entLine,
@@ -885,7 +886,7 @@ def test_initializeMassEnt():
 
     # call function to be tested
     simTypeActual = "res"
-    entrMassRaster, entrEnthRaster, reportAreaInfo = com1DFA.initializeMassEnt(
+    entrMassRaster, entrEnthRaster, entrDepthRaster, reportAreaInfo = com1DFA.initializeMassEnt(
         dem,
         simTypeActual,
         entLine,
@@ -2025,7 +2026,7 @@ def test_initializeFields():
     #    print("compute TA", fields["computeTA"])
     #    print("compute P", fields["computeP"])
 
-    assert len(fields) == 24
+    assert len(fields) == 25
     assert fields["computeTA"] is False
     assert fields["computeKE"] is False
     assert fields["computeP"]
@@ -2060,7 +2061,7 @@ def test_initializeFields():
     }
     # call function to be tested
     particles, fields = com1DFA.initializeFields(cfg, dem, particles, "")
-    assert len(fields) == 24
+    assert len(fields) == 25
     assert fields["computeTA"]
     assert fields["computeKE"]
     assert fields["computeP"] is False
@@ -2103,6 +2104,8 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "explicitFriction": 0,
         "timeDependentRelease": "False",
         "timeDependentReleaseScenarios": "",
+        "adaptSfcEntrainment": "0",
+        "entrainableDeposition": "False",
     }
     standardCfg["INPUT"] = {
         "entThThickness": "1.",
@@ -2179,6 +2182,8 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "explicitFriction": 0,
         "timeDependentRelease": "False",
         "timeDependentReleaseScenarios": "",
+        "adaptSfcEntrainment": "0",
+        "entrainableDeposition": "False",
     }
 
     testCfg["INPUT"] = {
@@ -2292,6 +2297,8 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "explicitFriction": 0,
         "timeDependentRelease": "False",
         "timeDependentReleaseScenarios": "",
+        "adaptSfcEntrainment": "0",
+        "entrainableDeposition": "False",
     }
     testCfg2["INPUT"] = {
         "entThThickness": "1.",
@@ -2371,6 +2378,8 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "rho": "200.0",
         "explicitFriction": 0,
         "timeDependentRelease": "True",
+        "adaptSfcEntrainment": "0",
+        "entrainableDeposition": "False",
     }
     standardCfg["INPUT"] = {
         "entThThickness": "1.",
@@ -2448,6 +2457,8 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "explicitFriction": 0,
         "timeDependentRelease": "True",
         "timeDependentReleaseScenarios": "release1PF",
+        "adaptSfcEntrainment": "0",
+        "entrainableDeposition": "False",
     }
 
     testCfg["INPUT"] = {
@@ -2482,6 +2493,18 @@ def test_prepareVarSimDict(tmp_path, caplog):
         for key in testCfg[section]:
             assert simDict[simName1]["cfgSim"][section][key] == testCfg[section][key]
 
+    standardCfg["GENERAL"]["entrainableDeposition"] = "True"
+    simDict = com1DFA.prepareVarSimDict(standardCfg, inputSimFiles, variationDict)
+
+    for key in testDict[simName1]:
+        #        print(simDict)
+        #        print(simDict[simName1][key])
+        assert simDict[simName1][key] == testDict[simName1][key]
+
+    for section in testCfg.sections():
+        for key in testCfg[section]:
+            assert simDict[simName1]["cfgSim"][section][key] == testCfg[section][key]
+
 
 def test_initializeSimulation(tmp_path):
     """test initializing a simulation"""
@@ -3163,6 +3186,7 @@ def test_adaptDEM():
         "adaptSfcStopped": 0,
         "adaptSfcDetrainment": 0,
         "adaptSfcEntrainment": 0,
+        "entrainableDeposition": "False",
     }
 
     header = {
@@ -3196,6 +3220,9 @@ def test_adaptDEM():
         "demAdapted": data,
         "sfcChangeTotal": np.zeros_like(data),
         "sfcChange": np.zeros_like(data),
+        "mStop": np.zeros_like(data),
+        "entrDepth": np.zeros_like(data),
+        "entrMassRaster": np.zeros_like(data),
     }
 
     dem = geoTrans.getNormalMesh(dem, num=cfg["GENERAL"].getfloat("methodMeshNormal"))
@@ -3218,6 +3245,7 @@ def test_adaptDEM():
         "adaptSfcStopped": 1,
         "adaptSfcDetrainment": 1,
         "adaptSfcEntrainment": 1,
+        "entrainableDeposition": "False",
     }
 
     # all rasters for depth changes are zero
@@ -3255,6 +3283,7 @@ def test_adaptDEM():
         "adaptSfcStopped": 1,
         "adaptSfcDetrainment": 0,
         "adaptSfcEntrainment": 1,
+        "entrainableDeposition": "False",
     }
 
     fieldsInput = fields.copy()
@@ -3272,6 +3301,7 @@ def test_adaptDEM():
         "adaptSfcStopped": 1,
         "adaptSfcDetrainment": 1,
         "adaptSfcEntrainment": 1,
+        "entrainableDeposition": "False",
     }
 
     fields["FTEnt"] -= 1
@@ -3319,6 +3349,65 @@ def test_adaptDEM():
     assert np.all(fieldsAdapted["sfcChange"] == fields["FTDet"] / NzNormed)
     assert np.all(fieldsAdapted["sfcChangeTotal"] == fields["FTDet"] / NzNormed)
 
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {
+        "methodMeshNormal": 1,
+        "adaptSfcStopped": 1,
+        "adaptSfcDetrainment": 0,
+        "adaptSfcEntrainment": 1,
+        "entrainableDeposition": "True",
+    }
+
+    fields["FTEnt"] = np.zeros_like(fields["FTDet"])
+    fields["FTStop"] = np.array(
+        [
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+            [1, 1, 1, 1, 1],
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+        ],
+        dtype=float,
+    )
+    fields["mStop"] = np.array(
+        [
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+            [-10, -10, -10, -10, -10],
+            [0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0],
+        ],
+        dtype=float,
+    )
+    fields["entrDepth"] = np.zeros_like(fields["FTDet"])
+    fields["entrMassRaster"] = np.zeros_like(fields["FTDet"])
+    fields["demNotErodableRaster"] = dem["rasterData"].copy()
+    fieldsInput = fields.copy()
+    demInput = dem.copy()
+    demAdapted, fieldsAdapted = com1DFA.adaptDEM(demInput, fieldsInput, cfg["GENERAL"])
+
+    assert np.all(
+        demAdapted["rasterData"]
+        == np.array(
+            [
+                [1.0, 2.0, 3.0, 4.0, 5.0],
+                [1.0, 2.0, 3.0, 4.0, 5.0],
+                [1.0, 2.0, 3.0, 4.0, 5.0] + 1 / NzNormed[2],
+                [1.0, 2.0, 3.0, 4.0, 5.0],
+                [1.0, 2.0, 3.0, 4.0, 5.0],
+            ]
+        )
+    )
+    assert np.any(demAdapted["Nx"] != dem["Nx"])
+    assert np.any(demAdapted["Ny"] != dem["Ny"])
+    assert np.all(demAdapted["Nz"] == dem["Nz"])
+    assert np.any(dem["areaRaster"] != demAdapted["areaRaster"])
+    assert np.all(fieldsAdapted["sfcChange"] == fields["FTDet"] / NzNormed)
+    assert np.all(fieldsAdapted["sfcChangeTotal"] == fields["FTDet"] / NzNormed)
+    assert np.all(fieldsAdapted["demNotErodableRaster"] == fields["demNotErodableRaster"])
+    assert np.all(fieldsAdapted["entrDepth"] == fields["FTStop"] / NzNormed)
+    assert np.all(fieldsAdapted["entrMassRaster"] == -fields["mStop"])
+
 
 def test_tSteps_output_behavior(tmp_path, caplog):
     """Test that tSteps controls which timesteps are exported correctly.