diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py
index 4ae908567..dbfa58313 100644
--- a/avaframe/com1DFA/com1DFA.py
+++ b/avaframe/com1DFA/com1DFA.py
@@ -1051,6 +1051,12 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     # create primary release area particles and fields
     releaseLine["header"] = dem["originalHeader"]
     inputSimLines["releaseLine"]["header"] = dem["originalHeader"]
+    # export release area raster to file
+    if cfg['EXPORTS'].getboolean('exportRasters'):
+        outDir = pathlib.Path(cfgGen['avalancheDir'], 'Outputs', 'internalRasters')
+        fU.makeADir(outDir)
+        IOf.writeResultToRaster(dem["originalHeader"], relRaster, (outDir / 'releaseRaster'), flip=True)
+        log.info('Release area raster derived from shp file saved to %s' % str(outDir / 'releaseRaster'))
     particles = initializeParticles(
         cfgGen,
         releaseLine,
@@ -1105,7 +1111,7 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
         inputSimLines["entLine"],
         reportAreaInfo,
         thresholdPointInPoly,
-        cfgGen,
+        cfg,
     )
 
     # check if entrainment and release overlap
@@ -1621,7 +1627,7 @@ def initializeSecRelease(inputSimLines, dem, relRaster, reportAreaInfo):
     return secondaryReleaseInfo, reportAreaInfo
 
 
-def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPointInPoly, cfgGen):
+def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPointInPoly, cfg):
     """Initialize mass for entrainment
 
     Parameters
@@ -1637,8 +1643,8 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin
     thresholdPointInPoly: float
         threshold val that decides if a point is in the polygon, on the line or
         very close but outside
-    cfgGen: config parser
-        General configuration
+    cfg: config parser
+        configuration
 
     Returns
     -------
@@ -1657,14 +1663,18 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin
         log.info("Entrainment area features: %s" % (entLine["Name"]))
         entLine = geoTrans.prepareArea(entLine, dem, thresholdPointInPoly, thList=entLine["thickness"])
         entrMassRaster = entLine["rasterData"]
+        if cfg['EXPORTS'].getboolean('exportRasters'):
+            outDir = pathlib.Path(cfg['GENERAL']['avalancheDir'], 'Outputs', 'internalRasters')
+            IOf.writeResultToRaster(dem["originalHeader"], entrMassRaster, (outDir / 'entrainmentRaster'), flip=True)
+            log.info('Release area raster derived from shp file saved to %s' % str(outDir / 'entrainmentRaster'))
         # ToDo: not used in samos but implemented
-        # tempRaster = cfgGen.getfloat('entTempRef') + (dem['rasterData'] - cfgGen.getfloat('entMinZ'))
-        # * cfgGen.getfloat('entTempGrad')
-        # entrEnthRaster = np.where(tempRaster < 0, tempRaster*cfgGen.getfloat('cpIce'),
-        #                           tempRaster*cfgGen.getfloat('cpWtr')/cfgGen.getfloat('hFusion'))
+        # 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,
-            cfgGen.getfloat("entTempRef") * cfgGen.getfloat("cpIce"),
+            cfg['GENERAL'].getfloat("entTempRef") * cfg['GENERAL'].getfloat("cpIce"),
             0,
         )
         reportAreaInfo["entrainment"] = "Yes"
@@ -1673,7 +1683,7 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin
         entrEnthRaster = np.zeros((nrows, ncols))
         reportAreaInfo["entrainment"] = "No"
 
-    entrMassRaster = entrMassRaster * cfgGen.getfloat("rhoEnt")
+    entrMassRaster = entrMassRaster * cfg['GENERAL'].getfloat("rhoEnt")
 
     return entrMassRaster, entrEnthRaster, reportAreaInfo
 
diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini
index 257a8b06e..366035dcc 100644
--- a/avaframe/com1DFA/com1DFACfg.ini
+++ b/avaframe/com1DFA/com1DFACfg.ini
@@ -473,4 +473,6 @@ unitpfv = ms-1
 # peak files and plots are exported, option to turn off exports when exportData is set to False
 # this affects export of peak files and also generation of peak file plots
 exportData = True
+# export release and optional entrainment raster files derived from shp files saved to Outputs/com1DFA/internalRasters
+exportRasters = False
 
diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py
index 36090bfe8..77be236d2 100644
--- a/avaframe/tests/test_com1DFA.py
+++ b/avaframe/tests/test_com1DFA.py
@@ -700,6 +700,7 @@ def test_initializeMassEnt():
         "cpIce": "2050.",
         "TIni": "-10.",
     }
+    cfg['EXPORTS'] = {'exportRasters': False}
 
     simTypeActual = "entres"
     dirName = pathlib.Path(__file__).parents[0]
@@ -727,7 +728,7 @@ def test_initializeMassEnt():
         entLine,
         reportAreaInfo,
         thresholdPointInPoly,
-        cfg["GENERAL"],
+        cfg,
     )
     testData = np.zeros((nrows, ncols))
     testData[0:11, 0:11] = 1.0 * 200.0
@@ -750,7 +751,7 @@ def test_initializeMassEnt():
         entLine,
         reportAreaInfo,
         thresholdPointInPoly,
-        cfg["GENERAL"],
+        cfg,
     )
 
     assert np.array_equal(entrMassRaster, np.zeros((nrows, ncols)))
@@ -2122,6 +2123,7 @@ def test_initializeSimulation(tmp_path):
         "TIni": "-10.",
         "ResistanceModel": "cRes",
     }
+    cfg['EXPORTS'] = {'exportRasters': "False"}
     # setup dem input
     demHeader = {}
     demHeader["xllcenter"] = 1.0
@@ -2291,6 +2293,7 @@ def test_initializeSimulation(tmp_path):
         "restitutionCoefficient": 1,
         "nIterDam": 1,
     }
+    cfg['EXPORTS'] = {'exportRasters': "False"}
     releaseLine = {
         "x": np.asarray([6.9, 8.5, 8.5, 6.9, 6.9]),
         "y": np.asarray([7.9, 7.9, 9.5, 9.5, 7.9]),