diff --git a/avaframe/com1DFA/DFAfunctionsCython.pxd b/avaframe/com1DFA/DFAfunctionsCython.pxd
index 7987a1417..22611f0f4 100644
--- a/avaframe/com1DFA/DFAfunctionsCython.pxd
+++ b/avaframe/com1DFA/DFAfunctionsCython.pxd
@@ -8,6 +8,8 @@
 
 cpdef (double, double) computeEntMassAndForce(double, double, double, double, double, double, double)
 
+cpdef (double, double, double) computeEntMassAndForceTjem(double, double, double, double, double, double, double)
+
 cpdef double computeDetMass(double, double, double, double)
 
 cpdef double computeResForce(double, double, double, double, int, int)
diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx
index 8738b850e..9b9c70030 100644
--- a/avaframe/com1DFA/DFAfunctionsCython.pyx
+++ b/avaframe/com1DFA/DFAfunctionsCython.pyx
@@ -27,7 +27,7 @@ import avaframe.in3Utils.geoTrans as geoTrans
 log = logging.getLogger(__name__)
 
 
-def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType):
+def computeForceC(cfg, particles, fields, dem, int frictType, int entrType, int resistanceType):
   """ compute forces acting on the particles (without the SPH component)
 
   Cython implementation implementation
@@ -44,6 +44,8 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
       dictionary with dem information
   frictType: int
     identifier for friction law to be used
+  entrType: int
+    identifier for entrainment law to be used
   resistanceType: int
     identifier for resistance model to be used
 
@@ -73,6 +75,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double BSamosAtMedium = cfg.getfloat('Bsamosatmedium')
   cdef double RSamosAtMedium = cfg.getfloat('Rsamosatmedium')
   cdef double entEroEnergy = cfg.getfloat('entEroEnergy')
+  cdef double entShearStrength = cfg.getfloat('entShearStrength')
   cdef double entShearResistance = cfg.getfloat('entShearResistance')
   cdef double entDefResistance = cfg.getfloat('entDefResistance')
   cdef double rho = cfg.getfloat('rho')
@@ -393,7 +396,10 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
       entrMassCell = entrMassRaster[indCellY, indCellX]
       if entrMassCell > 0.0:
         # compute entrained mass
-        dm, areaEntrPart = computeEntMassAndForce(dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt)
+        if entrType == 1:
+            dm, areaEntrPart = computeEntMassAndForce(dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt)
+        if entrType == 2:
+            dm, areaEntrPart, tau = computeEntMassAndForceTjem(dt, entrMassCell, areaPart, uMag, tau, entShearStrength, rhoEnt)
         # speed loss due to energy loss due to entrained mass
         dEnergyEntr = areaEntrPart * entShearResistance + dm * entDefResistance
 
@@ -542,6 +548,54 @@ cpdef (double, double) computeEntMassAndForce(double dt, double entrMassCell,
 
   return dm, areaEntrPart
 
+cpdef (double, double, double) computeEntMassAndForceTjem(double dt, double entrMassCell,
+                                              double areaPart, double uMag,
+                                              double tau, double entShearStrength,
+                                              double rhoEnt):
+  """ compute force component due to entrained mass
+
+  Parameters
+  ----------
+  dt: float
+    time step
+  entrMassCell : float
+      available mass for entrainement
+  areaPart : float
+      particle area
+  uMag : float
+      particle speed (velocity magnitude)
+  tau : float
+      bottom shear stress
+
+  Returns
+  -------
+  dm : float
+      entrained mass
+  areaEntrPart : float
+      Area for entrainment energy loss computation
+  entShearStrength: float
+    shear strength of the entrainable snow mass
+  rhoEnt: float
+    entrainement density
+  """
+  cdef double width, ABotSwiped
+  # compute entrained mass
+  cdef double dm = 0.0
+  cdef double areaEntrPart = 0.0
+
+  if entrMassCell > 0:
+
+      # erosion: erode according to shear strength of snow cover
+      # in contrast to MoTVoellmy tau_c here is not scaled with the flow density
+      # because com1DFA does not scale shear stress at the bottom of the avalanche with the flow density while MoTVoellmy does
+      tau_c = entShearStrength
+      dm = max(0, tau - tau_c) * areaPart * dt / uMag
+      areaEntrPart = areaPart
+      if dm > 0:
+        tau = min(tau, tau_c)
+
+  return dm, areaEntrPart, tau
+
 
 cpdef double computeDetMass(double dt, double detCell,
                                               double areaPart, double uMag):
diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py
index e928523d3..992203dfe 100644
--- a/avaframe/com1DFA/com1DFA.py
+++ b/avaframe/com1DFA/com1DFA.py
@@ -2189,6 +2189,17 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     frictType = frictModelsList.index(frictModel) + 1
     log.debug("Friction Model used: %s, %s" % (frictModelsList[frictType - 1], frictType))
 
+    # derive entrainment type
+    # turn entrainment model into integer
+    entrModelsList = [
+        "avaframe",
+        "tjem",
+    ]
+    entrModel = cfgGen["entrModel"].lower()
+    entrType = entrModelsList.index(entrModel) + 1
+    log.debug("Entrainment Model used: %s, %s" % (entrModelsList[entrType - 1], entrType))
+    print("Entrainment Model used: %s, %s" % (entrModelsList[entrType - 1], entrType))
+
     # turn resistance model into integer
     ResModel = cfgGen["ResistanceModel"].lower()
     ResModelsList = [
@@ -2301,6 +2312,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
             dem,
             tCPU,
             frictType,
+            entrType,
             resistanceType,
             inputSimLines["reportAreaInfo"],
         )
@@ -2723,6 +2735,7 @@ def computeEulerTimeStep(
     dem,
     tCPU,
     frictType,
+    entrType,
     resistanceType,
     reportAreaInfo,
 ):
@@ -2744,6 +2757,8 @@ def computeEulerTimeStep(
         computation time dictionary
     frictType: int
         indicator for chosen type of friction model
+    entrType: int
+        indicator for chosen type of entrainment model
     resistanceType: int
         identifier for chosen type of resistance model
 
@@ -2776,7 +2791,7 @@ def computeEulerTimeStep(
     startTime = time.time()
     # loop version of the compute force
     log.debug("Compute Force C")
-    particles, force, fields = DFAfunC.computeForceC(cfg, particles, fields, dem, frictType, resistanceType)
+    particles, force, fields = DFAfunC.computeForceC(cfg, particles, fields, dem, frictType, entrType, resistanceType)
     tCPUForce = time.time() - startTime
     tCPU["timeForce"] = tCPU["timeForce"] + tCPUForce
     # compute lateral force (SPH component of the calculation)
diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini
index 386282123..b5e9abeff 100644
--- a/avaframe/com1DFA/com1DFACfg.ini
+++ b/avaframe/com1DFA/com1DFACfg.ini
@@ -476,11 +476,23 @@ forestThMin = 0.6
 forestVMax = 40.
 forestThMax = 10.
 
+#++++++++++++ Entrainment Model
+### Available Models:
+#avaframe
+#tjem
+entrModel = avaframe
 
 #++++++++++++ Entrainment Erosion Energy
-# Used to determine speed loss via energy loss due to entrained mass
+# Used to determine speed loss via energy loss due to entrained mass when using Avaframe Entrainment Model
 # Select entEroEnergy = 0 to enable frontal plowing; entEroEnergy > 0 means basal erosion (without frontal plowing)
 entEroEnergy = 5000
+
+#++++++++++++ Entrainment Shear Strength
+# Used to determine the strength of the erodible snow cover when using TJEM Entrainment Model
+entShearStrength = 500
+
+#++++++++++++Entrainment Energy loss
+# Used to determine the energy loss because of entrainment
 entShearResistance = 0
 entDefResistance = 0