diff --git a/avaframe/ana4Stats/getStats.py b/avaframe/ana4Stats/getStats.py
index 71a79606d..4113dd87d 100644
--- a/avaframe/ana4Stats/getStats.py
+++ b/avaframe/ana4Stats/getStats.py
@@ -21,7 +21,16 @@
 log = logging.getLogger(__name__)
 
 
-def extractMaxValues(inputDir, avaDir, varPar, restrictType='', nameScenario='', parametersDict=''):
+def extractMaxValues(
+    inputDir,
+    avaDir,
+    varPar,
+    restrictType="",
+    nameScenario="",
+    parametersDict="",
+    layer="",
+    modName="com1DFA",
+):
     """ Extract max values of result parameters and save to dictionary
 
         - optionally restrict data of peak fields by defining which result parameter with restrictType,
@@ -51,12 +60,25 @@ def extractMaxValues(inputDir, avaDir, varPar, restrictType='', nameScenario='',
         """
 
     # filter simulation results using parametersDict
-    simNameList = cfgHandling.filterSims(avaDir, parametersDict)
+    simNameList = cfgHandling.filterSims(avaDir, parametersDict, modName=modName)
     # load dataFrame of all simulation configurations
-    simDF = cfgUtils.createConfigurationInfo(avaDir, standardCfg='')
+    simDF = cfgUtils.createConfigurationInfo(avaDir, standardCfg="", comModule=modName)
 
     # load peakFiles of all simulations and generate dataframe
     peakFilesDF = fU.makeSimDF(inputDir, avaDir=avaDir)
+    hasMultiLayer = any(peakFilesDF["layer"] != "")
+    if hasMultiLayer and layer == "":
+        message = (
+            "Multi-layer results detected but no layer specified. "
+            "Set 'layer' in getStatsCfg.ini (e.g. layer = L1)"
+        )
+        log.error(message)
+        raise ValueError(message)
+    elif not hasMultiLayer and layer != "":
+        message = "No multi-layer results detected but layer is specified. "
+        log.error(message)
+        raise ValueError(message)
+
     nSims = len(peakFilesDF['simName'])
     # initialize peakValues dictionary
     peakValues = {}
@@ -65,60 +87,79 @@ def extractMaxValues(inputDir, avaDir, varPar, restrictType='', nameScenario='',
 
     # Loop through result field files and compute statistical measures
     for m in range(nSims):
-        # filter simulations according to parametersDict
-        if peakFilesDF['simName'][m] in simNameList:
-
-            # Load data
-            fileName = peakFilesDF['files'][m]
-            simName = peakFilesDF['simName'][m]
-            dataFull = IOf.readRaster(fileName, noDataToNan=True)
-
-            # if restrictType, set result field values to nan if restrictType result field equals 0
-            if restrictType != '':
-                peakFilesSimName = peakFilesDF[peakFilesDF['simName'] == simName]
-                fileNameRestrict = peakFilesSimName[peakFilesSimName['resType'] == restrictType]['files'].values[0]
-                dataRestrict = IOf.readRaster(fileNameRestrict, noDataToNan=True)
-                data = np.where((dataRestrict['rasterData'] == 0.0), np.nan, dataFull['rasterData'])
-            else:
-                data = dataFull['rasterData']
-
-            # compute max, mean, min and standard deviation of result field
-            max = np.nanmax(data)
-            min = np.nanmin(data)
-            mean = np.nanmean(data)
-            std = np.nanstd(data)
-            statVals = {'max': max, 'min': min, 'mean': mean, 'std': std}
-            # add statistical measures
-            peakValues[simName].update({peakFilesDF['resType'][m]: statVals})
-
-            # fetch varPar value and nameScenario
-            varParVal = simDF[simDF['simName'] == simName][varPar].iloc[0]
-
-            # if varParVal is empty or nan check if thickness value
-            # if so - read from shp thickness value for each feature need to be found
-            if (varParVal == '' or pd.isna(varParVal)) and varPar in ['relTh', 'entTh', 'secondaryRelTh']:
-                # fetch id of thickness features
-                varParId = varPar + 'Id'
-                varParIds = str(simDF[simDF['simName'] == simName][varParId].iloc[0])
-                varParIdList = varParIds.split('|')
-                # create feature thickness parameter names
-                varParFirstName = varPar + varParIdList[0]
-                # chose value from first feature found!
-                varParVal = simDF[simDF['simName'] == simName][varParFirstName].iloc[0]
-                if len(varParIdList) > 1:
-                    log.warning('%s value - there are multiple features - %s value used '
-                                'from first feature only!' % (varPar, varPar))
-            # cadd varPar value to dictionary
-            peakValues[simName].update({'varPar': varParVal})
-            if nameScenario != '':
-                nameScenarioVal = simDF[simDF['simName'] == simName][nameScenario]
-                peakValues[simName].update({'scenario': nameScenarioVal.iloc[0]})
-                log.info('Simulation parameter %s= %s for resType: %s and name %s' %
-                        (varPar, str(varParVal), peakFilesDF['resType'][m], nameScenarioVal.iloc[0]))
-            else:
-                log.info('Simulation parameter %s= %s for resType: %s' %
-                        (varPar, str(varParVal), peakFilesDF['resType'][m]))
+        simName = peakFilesDF['simName'][m]
+        # skip simulations not in the filtered list
+        if simName not in simNameList:
+            peakValues.pop(simName, None)
+            continue
+
+        # Load data
+        fileName = peakFilesDF['files'][m]
+        dataFull = IOf.readRaster(fileName, noDataToNan=True)
+
+        # if restrictType, set result field values to nan if restrictType result field equals 0
+        if restrictType != "":
+            peakFilesSimName = peakFilesDF[peakFilesDF['simName'] == simName]
+            # mask for resultType
+            mask = peakFilesSimName["resType"] == restrictType
+            # mask for layer
+            if layer:
+                mask &= peakFilesSimName["layer"] == layer
+            fileNameRestrict1 = peakFilesSimName[mask]
+            if len(fileNameRestrict1) == 0:
+                message = "No %s file found for %s to restrict current result file %s with." % (
+                    restrictType,
+                    simName,
+                    fileName.name,
+                )
+                log.error(message)
+                raise ValueError(message)
+            fileNameRestrict = fileNameRestrict1["files"].to_list()[0]
+            dataRestrict = IOf.readRaster(fileNameRestrict, noDataToNan=True)
+            log.info("%s file restricted using %s." % (fileName.name, fileNameRestrict.name))
+            data = np.where((dataRestrict['rasterData'] == 0.0), np.nan, dataFull['rasterData'])
+        else:
+            data = dataFull['rasterData']
+
+        # compute max, mean, min and standard deviation of result field
+        max = np.nanmax(data)
+        min = np.nanmin(data)
+        mean = np.nanmean(data)
+        std = np.nanstd(data)
+        statVals = {'max': max, 'min': min, 'mean': mean, 'std': std}
+        # use layer-suffixed key for multi-layer results (e.g. ppr_l1)
+        resType = peakFilesDF['resType'][m]
+        if peakFilesDF['layer'][m]:
+            resType = "%s_%s" % (resType, peakFilesDF['layer'][m].lower())
+        peakValues[simName].update({resType: statVals})
+
+        # fetch varPar value and nameScenario
+        simRow = simDF[simDF['simName'] == simName]
+        varParVal = simRow[varPar].iloc[0]
+
+        # if varParVal is empty or nan check if thickness value
+        # if so - read from shp thickness value for each feature need to be found
+        if (varParVal == '' or pd.isna(varParVal)) and varPar in ['relTh', 'entTh', 'secondaryRelTh']:
+            # fetch id of thickness features
+            varParId = varPar + 'Id'
+            varParIds = str(simRow[varParId].iloc[0])
+            varParIdList = varParIds.split('|')
+            # create feature thickness parameter names
+            varParFirstName = varPar + varParIdList[0]
+            # chose value from first feature found!
+            varParVal = simRow[varParFirstName].iloc[0]
+            if len(varParIdList) > 1:
+                log.warning('%s value - there are multiple features - %s value used '
+                            'from first feature only!' % (varPar, varPar))
+        # add varPar value to dictionary
+        peakValues[simName].update({'varPar': varParVal})
+        if nameScenario != '':
+            nameScenarioVal = simRow[nameScenario]
+            peakValues[simName].update({'scenario': nameScenarioVal.iloc[0]})
+            log.info('Simulation parameter %s= %s for resType: %s and name %s' %
+                    (varPar, str(varParVal), peakFilesDF['resType'][m], nameScenarioVal.iloc[0]))
         else:
-            peakValues.pop(peakFilesDF['simName'][m], None)
+            log.info('Simulation parameter %s= %s for resType: %s' %
+                    (varPar, str(varParVal), peakFilesDF['resType'][m]))
 
-    return peakValues
+    return peakValues, hasMultiLayer
diff --git a/avaframe/ana4Stats/getStatsCfg.ini b/avaframe/ana4Stats/getStatsCfg.ini
index c5b9ff247..65e7083e9 100644
--- a/avaframe/ana4Stats/getStatsCfg.ini
+++ b/avaframe/ana4Stats/getStatsCfg.ini
@@ -20,6 +20,12 @@ varParUnit = m
 nameScenario = releaseScenario
 # restrict type for restricting the result field
 restrictType = pft
+# resultType 1 for plots
+resType1 = pft
+# resultType 2 for plots
+resType2 = pfv
+# if multilayer result files
+layer =
 
 
 ## Uncomment this section FILTER in your local copy of the ini file and add filter parameter and parameter values
diff --git a/avaframe/ana4Stats/probAna.py b/avaframe/ana4Stats/probAna.py
index eb7d5dc5f..18afc1dee 100644
--- a/avaframe/ana4Stats/probAna.py
+++ b/avaframe/ana4Stats/probAna.py
@@ -462,7 +462,9 @@ def probAnalysis(avaDir, cfg, modName, parametersDict="", inputDir="", probConf=
 
     # fetch all result files and filter simulations according to parametersDict
     if modName.lower() in ["com1dfa", "com5snowslide", "com6rockavalanche", "com8motpsa", "com9motvoellmy"]:
-        simNameList = cfgHandling.filterSims(avaDir, parametersDict, specDir=inputDir, simDF=simDFActual)
+        simNameList = cfgHandling.filterSims(
+            avaDir, parametersDict, specDir=inputDir, simDF=simDFActual, modName=modName
+        )
         filtering = True
     else:
         simNameList = []
@@ -490,6 +492,18 @@ def probAnalysis(avaDir, cfg, modName, parametersDict="", inputDir="", probConf=
 
         raise FileNotFoundError(message)
 
+    # check layer configuration for multi-layer results
+    peakVar = cfg["GENERAL"]["peakVar"]
+    layer = cfg["GENERAL"].get("layer", "")
+    hasMultiLayer = any(peakFilesDF["layer"] != "")
+    if hasMultiLayer and layer == "":
+        message = (
+            "Multi-layer results detected but no layer specified. "
+            "Set 'layer' in probAnaCfg.ini (e.g. layer = L1)"
+        )
+        log.error(message)
+        raise ValueError(message)
+
     # get header info from peak files - this should be the same for all peakFiles
     header = IOf.readRasterHeader(peakFilesDF["files"][0])
     refData = IOf.readRaster(peakFilesDF["files"][0])
@@ -505,8 +519,9 @@ def probAnalysis(avaDir, cfg, modName, parametersDict="", inputDir="", probConf=
     for m in range(len(peakFilesDF["names"])):
         # only take simulations that match filter criteria from parametersDict
         if (peakFilesDF["simName"][m] in simNameList) or filtering == False:
-            # Load peak field for desired peak field parameter
-            if peakFilesDF["resType"][m] == cfg["GENERAL"]["peakVar"]:
+            # Load peak field for desired peak field parameter, filtering by layer if set
+            layerMatch = (layer == "") or (peakFilesDF["layer"][m].lower() == layer.lower())
+            if peakFilesDF["resType"][m] == peakVar and layerMatch:
                 # Load data
                 fileName = peakFilesDF["files"][m]
                 dataLim = np.zeros((nRows, nCols))
@@ -534,7 +549,7 @@ def probAnalysis(avaDir, cfg, modName, parametersDict="", inputDir="", probConf=
                 )
                 contourDict[fileName.stem] = contourDictXY
 
-                log.info("File Name: %s , simulation parameter %s " % (fileName, cfg["GENERAL"]["peakVar"]))
+                log.info("File Name: %s , simulation parameter %s " % (fileName, peakVar))
 
                 # Check if peak values exceed desired threshold
                 dataLim[data > float(cfg["GENERAL"]["peakLim"])] = 1.0
@@ -543,19 +558,21 @@ def probAnalysis(avaDir, cfg, modName, parametersDict="", inputDir="", probConf=
 
     # Create probability map ranging from 0-1
     probMap = probSum / count
-    unit = pU.cfgPlotUtils["unit%s" % cfg["GENERAL"]["peakVar"]]
+    unit = pU.cfgPlotUtils["unit%s" % peakVar]
+    layerStr = "_%s" % layer if layer else ""
     log.info(
-        "probability analysis performed for peak parameter: %s and a peak value "
-        "threshold of: %s %s" % (cfg["GENERAL"]["peakVar"], cfg["GENERAL"]["peakLim"], unit)
+        "probability analysis performed for peak parameter: %s, layer: %s and a peak value "
+        "threshold of: %s %s" % (peakVar, layer or "single-layer", cfg["GENERAL"]["peakLim"], unit)
     )
     log.info("%s peak fields added to analysis" % count)
 
     # Save to raster file
     avaName = avaDir.name
-    outFileName = "%s_prob_%s_%s_lim%s" % (
+    outFileName = "%s_prob_%s_%s%s_lim%s" % (
         avaName,
         probConf,
-        cfg["GENERAL"]["peakVar"],
+        peakVar,
+        layerStr,
         cfg["GENERAL"]["peakLim"],
     )
     outFile = outDir / outFileName
diff --git a/avaframe/ana4Stats/probAnaCfg.ini b/avaframe/ana4Stats/probAnaCfg.ini
index 9f2d67f25..1859b48b4 100644
--- a/avaframe/ana4Stats/probAnaCfg.ini
+++ b/avaframe/ana4Stats/probAnaCfg.ini
@@ -6,6 +6,8 @@
 [GENERAL]
 # desired result parameter
 peakVar = ppr
+# layer for multi-layer results (e.g. L1, L2), leave empty for single-layer
+layer =
 # unit of desired result parameter
 unit = kPa
 # threshold for probability computations
diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py
index 472c0fda9..cc7e109a9 100644
--- a/avaframe/com1DFA/com1DFA.py
+++ b/avaframe/com1DFA/com1DFA.py
@@ -90,6 +90,13 @@ def com1DFAPreprocess(cfgMain, cfgInfo, module=com1DFA):
 
     avalancheDir = cfgMain["MAIN"]["avalancheDir"]
 
+    # batch mode: cfgInfo is a directory containing multiple .ini files
+    if isinstance(cfgInfo, pathlib.Path) and cfgInfo.is_dir():
+        simDict, inputSimFiles, simDFExisting, outDir = com1DFATools.createSimDictFromCfgs(
+            cfgMain, cfgInfo, module=module
+        )
+        return simDict, outDir, inputSimFiles, simDFExisting
+
     # read initial configuration
     if isinstance(cfgInfo, configparser.ConfigParser):
         cfgStart = cfgInfo
@@ -136,13 +143,7 @@ def com1DFAMain(cfgMain, cfgInfo=""):
 
     avalancheDir = cfgMain["MAIN"]["avalancheDir"]
 
-    # Route based on cfgInfo type: directory Path = batch mode, otherwise single config mode
-    if isinstance(cfgInfo, pathlib.Path) and cfgInfo.is_dir():
-        # Batch mode - cfgInfo is directory path (from getModuleConfig with batchCfgDir)
-        simDict, inputSimFiles, simDFExisting, outDir = com1DFATools.createSimDictFromCfgs(cfgMain, cfgInfo)
-    else:
-        # Single config mode - cfgInfo is ConfigParser, file path string, or empty string
-        simDict, outDir, inputSimFiles, simDFExisting = com1DFAPreprocess(cfgMain, cfgInfo)
+    simDict, outDir, inputSimFiles, simDFExisting = com1DFAPreprocess(cfgMain, cfgInfo)
 
     log.info("The following simulations will be performed")
     for key in simDict:
diff --git a/avaframe/com8MoTPSA/com8MoTPSA.py b/avaframe/com8MoTPSA/com8MoTPSA.py
index 5cdad5418..033e7d4ec 100644
--- a/avaframe/com8MoTPSA/com8MoTPSA.py
+++ b/avaframe/com8MoTPSA/com8MoTPSA.py
@@ -35,7 +35,7 @@ def com8MoTPSAMain(cfgMain, cfgInfo=None):
     """
     # Get all necessary information from the configuration files
     currentModule = sys.modules[__name__]
-    simDict, inputSimFiles = mT.MoTGenerateConfigs(cfgMain, cfgInfo, currentModule)
+    simDict, _, inputSimFiles, _ = com1DFA.com1DFAPreprocess(cfgMain, cfgInfo, module=currentModule)
 
     # convert DEM from nan to 0 values
     # TODO: suggest MoT-PSA to handle nan values
diff --git a/avaframe/com8MoTPSA/runAna4ProbAnaCom8MoTPSA.py b/avaframe/com8MoTPSA/runAna4ProbAnaCom8MoTPSA.py
index bc8c874ad..0d482bb7e 100644
--- a/avaframe/com8MoTPSA/runAna4ProbAnaCom8MoTPSA.py
+++ b/avaframe/com8MoTPSA/runAna4ProbAnaCom8MoTPSA.py
@@ -110,7 +110,8 @@ def runProbAna(avalancheDir=''):
                     }
         oP.plotContours(contourDict,
                         cfgProb['GENERAL']['peakVar'],
-                        cfgProb['GENERAL']['peakLim'], pathDict, addLegend=False)
+                        cfgProb['GENERAL']['peakLim'], pathDict, addLegend=False,
+                        layer=cfgProb['GENERAL'].get('layer', ''))
 
     # plot probability maps
     sP.plotProbMap(avalancheDir, inputDir, cfgProb, demPlot=True)
diff --git a/avaframe/com9MoTVoellmy/com9MoTVoellmy.py b/avaframe/com9MoTVoellmy/com9MoTVoellmy.py
index 5068a37ec..b6f621174 100644
--- a/avaframe/com9MoTVoellmy/com9MoTVoellmy.py
+++ b/avaframe/com9MoTVoellmy/com9MoTVoellmy.py
@@ -54,7 +54,7 @@ def com9MoTVoellmyMain(cfgMain, cfgInfo=None):
 
     # Get all necessary information from the configuration files
     currentModule = sys.modules[__name__]  # As if you would to import com9MoTVoellmy
-    simDict, inputSimFiles = mT.MoTGenerateConfigs(cfgMain, cfgInfo, currentModule)
+    simDict, _, inputSimFiles, _ = com1DFA.com1DFAPreprocess(cfgMain, cfgInfo, module=currentModule)
 
     # convert DEM from nan to 0 values
     # TODO: suggest MoT-PSA to handle nan values
diff --git a/avaframe/in3Utils/MoTUtils.py b/avaframe/in3Utils/MoTUtils.py
index 89b7278c8..d04d8126d 100644
--- a/avaframe/in3Utils/MoTUtils.py
+++ b/avaframe/in3Utils/MoTUtils.py
@@ -7,7 +7,6 @@
 
 import numpy as np
 
-from avaframe.com1DFA import com1DFA as com1DFA
 from avaframe.in2Trans import rasterUtils as rU
 
 log = logging.getLogger(__name__)
@@ -124,38 +123,6 @@ def runAndCheckMoT(command):
                 log.info(line)
 
 
-def MoTGenerateConfigs(cfgMain, cfgInfo, currentModule):
-    """
-    Creates configuration objects for com8MoTPSA and com9MoTVoellmy.
-
-    Parameters
-    ------------
-    cfgMain: configparser object
-        main configuration of AvaFrame
-    cfgInfo: str or pathlib Path or configparser object
-        path to configuration file if overwrite is desired - optional
-        if not local (if available) or default configuration will be loaded
-        if cfgInfo is a configparser object take this as initial config
-    currentModule: module object
-        is being passed to cfgUtils.writeCfgFile to create the correct the cfg
-
-
-    Returns
-    --------
-    simDict: dict
-        dictionary with one key per simulation to perform including its config object
-    inputSimFiles: dict
-        dictionary with input files info
-    """
-
-    # preprocessing to create configuration objects for all simulations to run
-    simDict, outDir, inputSimFiles, simDFExisting = com1DFA.com1DFAPreprocess(
-        cfgMain, cfgInfo, module=currentModule
-    )
-
-    return simDict, inputSimFiles
-
-
 def copyMoTFiles(workDir, outputDir, searchString, replaceString):
     """
     Copy and rename MoT result files from work directory to output directory.
diff --git a/avaframe/in3Utils/cfgHandling.py b/avaframe/in3Utils/cfgHandling.py
index 730b95dd3..61cfc1aaa 100644
--- a/avaframe/in3Utils/cfgHandling.py
+++ b/avaframe/in3Utils/cfgHandling.py
@@ -93,7 +93,7 @@ def addInfoToSimName(avalancheDir, csvString=""):
     return simDF[vars]
 
 
-def filterSims(avalancheDir, parametersDict, specDir="", simDF=""):
+def filterSims(avalancheDir, parametersDict, specDir="", simDF="", modName='com1DFA'):
     """Filter simulations using a list of parameters and a pandas dataFrame of simulation configurations
     if ~ is used as a prefix for a parameter - it is filtered according to values that do NOT match the value
     provided with the ~Parameter
@@ -108,6 +108,8 @@ def filterSims(avalancheDir, parametersDict, specDir="", simDF=""):
         path to a directory where simulation configuration files can be found - optional
     simDF: pandas DataFrame
         optional - if simDF already available
+    modName: str
+        name of computational module used to perform simulations to be filtered
 
     Returns
     --------
@@ -118,7 +120,7 @@ def filterSims(avalancheDir, parametersDict, specDir="", simDF=""):
     if isinstance(simDF, pd.DataFrame) is False:
         # load dataFrame for all configurations
         simDF = cfgUtils.createConfigurationInfo(
-            avalancheDir, standardCfg="", writeCSV=False, specDir=specDir
+            avalancheDir, standardCfg="", writeCSV=False, specDir=specDir, comModule=modName
         )
 
     # filter simulations all conditions in the parametersDict have to be met
diff --git a/avaframe/in3Utils/fileHandlerUtils.py b/avaframe/in3Utils/fileHandlerUtils.py
index f2ea2fcb7..d15a35d8a 100644
--- a/avaframe/in3Utils/fileHandlerUtils.py
+++ b/avaframe/in3Utils/fileHandlerUtils.py
@@ -589,6 +589,7 @@ def makeSimDF(inputDir, avaDir="", simID="simID"):
         "cellSize": [],
         simID: [],
         "timeStep": [],
+        "layer": [],
     }
 
     # Set name of avalanche if avaDir is given
@@ -614,6 +615,7 @@ def makeSimDF(inputDir, avaDir="", simID="simID"):
         data["modelType"].append(simNameParts["modelType"])
         data["resType"].append(simNameParts["resType"] if simNameParts["resType"] else "")
         data["timeStep"].append(simNameParts["timeStep"] if simNameParts["timeStep"] else "")
+        data["layer"].append(simNameParts["layer"] if simNameParts["layer"] else "")
 
         # Reconstruct simName (without resType and timeStep)
         # Preserve _AF_ separator if present in original name
diff --git a/avaframe/out3Plot/outQuickPlot.py b/avaframe/out3Plot/outQuickPlot.py
index 1b2e7aac6..9bfa30925 100644
--- a/avaframe/out3Plot/outQuickPlot.py
+++ b/avaframe/out3Plot/outQuickPlot.py
@@ -547,7 +547,7 @@ def generateOnePlot(dataDict, outDir, cfg, plotDict):
     return plotDict
 
 
-def plotContours(contourDict, resType, thresholdValue, pathDict, addLegend=True):
+def plotContours(contourDict, resType, thresholdValue, pathDict, addLegend=True, layer=""):
     """ plot contour lines of all transformed fields
 
         Parameters
@@ -562,14 +562,17 @@ def plotContours(contourDict, resType, thresholdValue, pathDict, addLegend=True)
             dictionary with info on project name, ...
         addLegend: bool
             if True add legend to plot
+        layer: str
+            layer identifier for multi-layer results (e.g. L1, L2), empty for single-layer
     """
 
     unit = pU.cfgPlotUtils['unit' + resType]
+    resTypeLabel = '%s (%s)' % (resType, layer) if layer else resType
     fig = plt.figure(figsize=(pU.figW*2, pU.figH))
 
     # show flow path
     ax1 = fig.add_subplot(121)
-    ax1.set_title('%s %s %s contour lines' % (resType, thresholdValue, unit))
+    ax1.set_title('%s %s %s contour lines' % (resTypeLabel, thresholdValue, unit))
     ax1.set_ylabel('x [m]')
     ax1.set_xlabel('y [m]')
     pU.putAvaNameOnPlot(ax1, pathDict['avaDir'])
diff --git a/avaframe/out3Plot/statsPlots.py b/avaframe/out3Plot/statsPlots.py
index 42c14541d..98fa163bf 100644
--- a/avaframe/out3Plot/statsPlots.py
+++ b/avaframe/out3Plot/statsPlots.py
@@ -10,6 +10,7 @@
 import seaborn as sns
 import pandas as pd
 import pathlib
+import re
 
 # local imports
 import avaframe.out3Plot.plotUtils as pU
@@ -26,26 +27,30 @@
 log = logging.getLogger(__name__)
 
 
-def plotValuesScatter(peakValues, resType1, resType2, cfg, avalancheDir, statsMeasure='max', flagShow=False):
-    """ Produce scatter plot of statistical measures (eg. max values) (resType1 und resType2),
-        for one set of simulations or multiple
+def plotValuesScatter(
+    peakValues, resType1, resType2, cfg, avalancheDir, statsMeasure="max", flagShow=False, layer=""
+):
+    """Produce scatter plot of statistical measures (eg. max values) (resType1 und resType2),
+    for one set of simulations or multiple
 
-        Parameters
-        -----------
-        peakDictList: dict
-            peakValues dictionary that contain max values of peak parameters and parameter variation info
-        resType1: str
-            result parameter 1, 'ppr', 'pft', 'pfv'
-        resType2: str
-            result parameter 1, 'ppr', 'pft', 'pfv'
-        cfg: dict
-            configuration, for now contains output location and varPar: parameter that is varied
-            to perfom a set of simulations
-        statsMeasure: str
-            statistical measure for plotting, options: max, mean, min, std
-        flagShow: bool
-            if True show plot
-        """
+    Parameters
+    -----------
+    peakDictList: dict
+        peakValues dictionary that contain max values of peak parameters and parameter variation info
+    resType1: str
+        result parameter 1, 'ppr', 'pft', 'pfv'
+    resType2: str
+        result parameter 1, 'ppr', 'pft', 'pfv'
+    cfg: dict
+        configuration, for now contains output location and varPar: parameter that is varied
+        to perfom a set of simulations
+    statsMeasure: str
+        statistical measure for plotting, options: max, mean, min, std
+    flagShow: bool
+        if True show plot
+    layer: str
+        layer info for multilayer results
+    """
 
     varPar = cfg['varPar']
     # extract values from dictionaries
@@ -60,10 +65,17 @@ def plotValuesScatter(peakValues, resType1, resType2, cfg, avalancheDir, statsMe
     log.info('Number of simulations is: %d' % (len(varVal)))
 
     # Get name and units for resTypes and parameter variation to annotate plots
-    name1 = pU.cfgPlotUtils['name%s' % resType1]
-    name2 = pU.cfgPlotUtils['name%s' % resType2]
-    unit1 = pU.cfgPlotUtils['unit%s' % resType1]
-    unit2 = pU.cfgPlotUtils['unit%s' % resType2]
+    baseResType1 = re.sub(r"_l\d+$", "", resType1)
+    baseResType2 = re.sub(r"_l\d+$", "", resType2)
+    # if multilayer add layer info to label names
+    if layer != "":
+        suffix = " (%s)" % layer
+    else:
+        suffix = ""
+    name1 = pU.cfgPlotUtils["name%s" % baseResType1] + suffix
+    name2 = pU.cfgPlotUtils["name%s" % baseResType2] + suffix
+    unit1 = pU.cfgPlotUtils["unit%s" % baseResType1]
+    unit2 = pU.cfgPlotUtils["unit%s" % baseResType2]
     nameVar = cfg['varParName']
     unitVar = cfg['varParUnit']
 
@@ -95,25 +107,36 @@ def plotValuesScatter(peakValues, resType1, resType2, cfg, avalancheDir, statsMe
     plotPath = pU.saveAndOrPlot({'pathResult': outDir}, plotName, fig)
 
 
-def plotValuesScatterHist(peakValues, resType1, resType2, cfg, avalancheDir,
-                          statsMeasure='max', flagShow=False, flagHue=False):
-    """ Produce scatter and marginal kde plot of max values, for one set of simulations or multiple
+def plotValuesScatterHist(
+    peakValues,
+    resType1,
+    resType2,
+    cfg,
+    avalancheDir,
+    statsMeasure="max",
+    flagShow=False,
+    flagHue=False,
+    layer="",
+):
+    """Produce scatter and marginal kde plot of max values, for one set of simulations or multiple
 
-        Parameters
-        -----------
-        peakValues: dict
-            peakValues dictionary that contain max values of peak parameters and parameter variation info
-        resType1: str
-            result parameter 1, 'ppr', 'pft', 'pfv'
-        resType2: str
-            result parameter 1, 'ppr', 'pft', 'pfv'
-        cfg: dict
-            configuration, for now contains output location and varPar: parameter that is varied
-            to perfom a set of simulations
-        statsMeasure: str
-            statistical measure for plotting, options: max, mean, min, std
-        flagShow: bool
-            if True show plot
+    Parameters
+    -----------
+    peakValues: dict
+        peakValues dictionary that contain max values of peak parameters and parameter variation info
+    resType1: str
+        result parameter 1, 'ppr', 'pft', 'pfv'
+    resType2: str
+        result parameter 1, 'ppr', 'pft', 'pfv'
+    cfg: dict
+        configuration, for now contains output location and varPar: parameter that is varied
+        to perfom a set of simulations
+    statsMeasure: str
+        statistical measure for plotting, options: max, mean, min, std
+    flagShow: bool
+        if True show plot
+    layer: str
+        layer info for multilayer results
 
     """
 
@@ -123,7 +146,7 @@ def plotValuesScatterHist(peakValues, resType1, resType2, cfg, avalancheDir,
     values1 = []
     values2 = []
     scenario = []
-
+    # if multilayer add layer info to label names
     for key in peakValues:
         values1.append(peakValues[key][resType1][statsMeasure])
         values2.append(peakValues[key][resType2][statsMeasure])
@@ -134,10 +157,16 @@ def plotValuesScatterHist(peakValues, resType1, resType2, cfg, avalancheDir,
     log.info('Number of simulations is: %d' % (len(varVal)))
 
     # Get name and units for resTypes and parameter variation to annotate plots
-    name1 = pU.cfgPlotUtils['name%s' % resType1]
-    name2 = pU.cfgPlotUtils['name%s' % resType2]
-    unit1 = pU.cfgPlotUtils['unit%s' % resType1]
-    unit2 = pU.cfgPlotUtils['unit%s' % resType2]
+    baseResType1 = re.sub(r"_l\d+$", "", resType1)
+    baseResType2 = re.sub(r"_l\d+$", "", resType2)
+    if layer != "":
+        suffix = " (%s)" % layer
+    else:
+        suffix = ""
+    name1 = pU.cfgPlotUtils["name%s" % baseResType1] + suffix
+    name2 = pU.cfgPlotUtils["name%s" % baseResType2] + suffix
+    unit1 = pU.cfgPlotUtils["unit%s" % baseResType1]
+    unit2 = pU.cfgPlotUtils["unit%s" % baseResType2]
     nameVar = cfg['varParName']
     unitVar = cfg['varParUnit']
     varValV = np.array(varVal)
@@ -309,9 +338,11 @@ def plotProbMap(avaDir, inDir, cfgFull, demPlot=False):
         # create figure and add title
         fig = plt.figure(figsize=(pU.figW*2, pU.figH))
 
+        layer = cfgFull['GENERAL'].get('layer', '')
+        peakVarLabel = '%s (%s)' % (cfgFull['GENERAL']['peakVar'], layer) if layer else cfgFull['GENERAL']['peakVar']
         fullTitle = '%s %s based on %s $>$ %s %s' % (avaName,
                                                      cfg['name'],
-                                                     cfgFull['GENERAL']['peakVar'],
+                                                     peakVarLabel,
                                                      cfgFull['GENERAL']['peakLim'],
                                                      cfgFull['GENERAL']['unit'])
         suptitle = fig.suptitle(fullTitle, fontsize=14, color='0.5')
diff --git a/avaframe/runAna4ProbAna.py b/avaframe/runAna4ProbAna.py
index a75672158..3f567f7be 100644
--- a/avaframe/runAna4ProbAna.py
+++ b/avaframe/runAna4ProbAna.py
@@ -106,7 +106,8 @@ def runProbAna(avalancheDir=''):
                     }
         oP.plotContours(contourDict,
                         cfgProb['GENERAL']['peakVar'],
-                        cfgProb['GENERAL']['peakLim'], pathDict, addLegend=False)
+                        cfgProb['GENERAL']['peakLim'], pathDict, addLegend=False,
+                        layer=cfgProb['GENERAL'].get('layer', ''))
 
     # plot probability maps
     sP.plotProbMap(avalancheDir, inputDir, cfgProb, demPlot=True)
diff --git a/avaframe/runProbAnalysisOnly.py b/avaframe/runProbAnalysisOnly.py
index 5c2da082a..557227234 100644
--- a/avaframe/runProbAnalysisOnly.py
+++ b/avaframe/runProbAnalysisOnly.py
@@ -79,7 +79,8 @@ def runProbAna(avalancheDir="", modName=""):
     )
     pathDict = {"pathResult": str(inputDir / "plots"), "avaDir": str(avalancheDir), "plotScenario": outName}
     oP.plotContours(
-        contourDict, cfgProb["GENERAL"]["peakVar"], cfgProb["GENERAL"]["peakLim"], pathDict, addLegend=False
+        contourDict, cfgProb["GENERAL"]["peakVar"], cfgProb["GENERAL"]["peakLim"], pathDict, addLegend=False,
+        layer=cfgProb["GENERAL"].get("layer", "")
     )
 
     # plot probability maps
diff --git a/avaframe/runScripts/runPlotProfile.py b/avaframe/runScripts/runPlotProfile.py
index ff4cad96f..c09642ab0 100644
--- a/avaframe/runScripts/runPlotProfile.py
+++ b/avaframe/runScripts/runPlotProfile.py
@@ -5,6 +5,7 @@
 # Load modules
 # importing general python modules
 import pathlib
+import re
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.cm import ScalarMappable
@@ -18,10 +19,11 @@
 from avaframe.out3Plot import plotUtils as pU
 import avaframe.in2Trans.rasterUtils as IOf
 
-
 ################USER Input#############
 varPar = "relTh0"
 resType = "pfv"
+layer = ""
+modName = "com1DFA"
 ############################################################
 
 # Load avalanche directory from general configuration file
@@ -37,13 +39,29 @@
 log.info("Current avalanche: %s", avalancheDir)
 
 # load dataFrame for all configurations
-simDF = cfgUtils.createConfigurationInfo(avalancheDir)
+simDF = cfgUtils.createConfigurationInfo(avalancheDir, comModule=modName)
 
 # create data frame that lists all available simulations and path to their result type result files
-inputsDF, resTypeList = fU.makeSimFromResDF(avalancheDir, "com1DFA")
+inputsDF, resTypeList = fU.makeSimFromResDF(avalancheDir, comModule=modName)
 # merge  parameters as columns to dataDF for matching simNames
 dataDF = inputsDF.merge(simDF, left_on="simName", right_on="simName")
 
+# check for multilayer results
+# fetch all result files
+peakFilesDF = fU.makeSimDF(pathlib.Path(avalancheDir, "Outputs", modName, "peakFiles"), avaDir=avalancheDir)
+hasMultiLayer = any(peakFilesDF["layer"] != "")
+if hasMultiLayer and layer == "":
+    message = (
+        "Multi-layer results detected but no layer specified. "
+        "Set 'layer' in the beginning of the script (e.g. layer = L1)"
+    )
+    log.error(message)
+    raise ValueError(message)
+elif hasMultiLayer and layer != "":
+    resTypeFile = resType + "_" + layer.lower()
+else:
+    resTypeFile = resType
+
 # fetch unique values of varPar
 values = sorted(set(dataDF[varPar].to_list()))
 log.info('%s values used for colorcoding are: %s' % (varPar, values))
@@ -51,20 +69,20 @@
 maxVal = np.nanmax(values)
 cmapSCVals = np.linspace(0, 1, len(values))
 # create colormap and setup ticks and itemsList
-unit = pU.cfgPlotUtils['unit' + resType]
+unit = pU.cfgPlotUtils["unit" + resType]
 cmapSC, colorSC, ticksSC, normSC, unitSC, itemsList, displayColorBar = pU.getColors4Scatter(values, len(values),
                                                                                                 unit)
 # create figure
 fig = plt.figure(figsize=(pU.figW * 2, pU.figH * 1.5))
 ax1 = fig.add_subplot(211)
 ax2 = fig.add_subplot(212)
-testField = IOf.readRaster(dataDF['pfv'].iloc[0])
+testField = IOf.readRaster(dataDF[resTypeFile].iloc[0])
 nrowsHalf = int(testField['rasterData'].shape[0] * 0.5)
 log.info('Dimension of raster: %d, %d' % (testField['rasterData'].shape[0], testField['rasterData'].shape[1]))
 log.info('Profile plotted for row %d' % nrowsHalf)
 for index, row in dataDF.iterrows():
     # read field
-    field = IOf.readRaster(row[resType])
+    field = IOf.readRaster(row[resTypeFile])
     fieldData = field['rasterData']
     cmapVal = cmapSCVals[values.index(row[varPar])]
     ax1.plot(fieldData[nrowsHalf, :],  c=cmapSC(cmapVal))
@@ -80,4 +98,3 @@
 ax2.hlines(nrowsHalf, 0, testField['rasterData'].shape[1], "white", label="profile location")
 ax2.legend()
 pU.saveAndOrPlot({"pathResult": outDir}, ('profile_%s_%s' % (resType, varPar)), fig)
-
diff --git a/avaframe/runScripts/runProbAna.py b/avaframe/runScripts/runProbAna.py
index 7998549a3..716540cfa 100644
--- a/avaframe/runScripts/runProbAna.py
+++ b/avaframe/runScripts/runProbAna.py
@@ -91,7 +91,8 @@
                     }
         oP.plotContours(contourDict,
                         cfgProb['GENERAL']['peakVar'],
-                        cfgProb['GENERAL']['peakLim'], pathDict)
+                        cfgProb['GENERAL']['peakLim'], pathDict,
+                        layer=cfgProb['GENERAL'].get('layer', ''))
 
     # plot probability maps
     sP.plotProbMap(avaDir, inputDir, cfgProb, demPlot=True)
diff --git a/avaframe/runScripts/runStatsExample.py b/avaframe/runScripts/runStatsExample.py
index c3137468c..f13bf2af0 100644
--- a/avaframe/runScripts/runStatsExample.py
+++ b/avaframe/runScripts/runStatsExample.py
@@ -92,9 +92,14 @@
 parametersDict = fU.getFilterDict(cfgStats, 'FILTER')
 
 # get statisical measure of simulations
-peakValues = getStats.extractMaxValues(inputDir, avaDir, cfgStats['GENERAL']['varPar'],
-                      restrictType=cfgStats['GENERAL']['restrictType'], nameScenario=cfgStats['GENERAL']['nameScenario'],
-                      parametersDict=parametersDict)
+peakValues, _ = getStats.extractMaxValues(
+    inputDir,
+    avaDir,
+    cfgStats["GENERAL"]["varPar"],
+    restrictType=cfgStats["GENERAL"]["restrictType"],
+    nameScenario=cfgStats["GENERAL"]["nameScenario"],
+    parametersDict=parametersDict,
+)
 
 # log to screen
 for key in peakValues:
diff --git a/avaframe/runScripts/runStatsExampleSimple.py b/avaframe/runScripts/runStatsExampleSimple.py
new file mode 100644
index 000000000..95a4558e0
--- /dev/null
+++ b/avaframe/runScripts/runStatsExampleSimple.py
@@ -0,0 +1,92 @@
+"""
+Run script for computing statistics of results including the runout derived from aimec
+"""
+
+# Load modules
+import pathlib
+
+
+# Local imports
+from avaframe.out3Plot import statsPlots as sPlot
+from avaframe.ana4Stats import getStats
+from avaframe.in3Utils import fileHandlerUtils as fU
+from avaframe.in3Utils import cfgUtils
+from avaframe.in3Utils import logUtils
+
+
+# log file name; leave empty to use default runLog.log
+logName = "runGetStats"
+
+# ++++++Set comModule name+++++++++
+modName = "com1DFA"
+
+# Load general configuration filee
+cfgMain = cfgUtils.getGeneralConfig()
+flagShow = cfgMain["FLAGS"].getboolean("showPlot")
+
+avaDir = cfgMain["MAIN"]["avalancheDir"]
+cfgStats = cfgUtils.getModuleConfig(getStats, avaDir)
+cfg = cfgStats["GENERAL"]
+
+# set output directory, first ava in list
+avaDir = pathlib.Path(avaDir)
+outDir = avaDir / "Outputs" / "ana4Stats"
+cfgStats["GENERAL"]["outDir"] = str(outDir)
+# Specify where you want the results to be stored
+fU.makeADir(outDir)
+
+# Start logging
+log = logUtils.initiateLogger(outDir, logName)
+
+# ----- determine max values of peak fields
+# set directory of peak files
+inputDir = avaDir / "Outputs" / modName / "peakFiles"
+
+# provide optional filter criteria for simulations
+parametersDict = fU.getFilterDict(cfgStats, "FILTER")
+
+# get statistical measure of simulations
+peakValues, hasMultiLayer = getStats.extractMaxValues(
+    inputDir,
+    avaDir,
+    cfgStats["GENERAL"]["varPar"],
+    restrictType=cfgStats["GENERAL"]["restrictType"],
+    nameScenario=cfgStats["GENERAL"]["nameScenario"],
+    parametersDict=parametersDict,
+    layer=cfgStats["GENERAL"].get("layer", ""),
+    modName=modName,
+)
+
+
+# ++++++++++++++ Plot max values +++++++++++++++++
+# create resType + layer for fetching the results if multilayer
+if hasMultiLayer:
+    resType1 = cfgStats["GENERAL"]["resType1"] + "_" + cfgStats["GENERAL"]["layer"].lower()
+    resType2 = cfgStats["GENERAL"]["resType2"] + "_" + cfgStats["GENERAL"]["layer"].lower()
+else:
+    resType1 = cfgStats["GENERAL"]["resType1"]
+    resType2 = cfgStats["GENERAL"]["resType2"]
+
+sPlot.plotValuesScatter(
+    peakValues,
+    resType1,
+    resType2,
+    cfgStats["GENERAL"],
+    avaDir,
+    statsMeasure="max",
+    flagShow=flagShow,
+    layer=cfgStats["GENERAL"]["layer"],
+)
+sPlot.plotValuesScatterHist(
+    peakValues,
+    resType1,
+    resType2,
+    cfgStats["GENERAL"],
+    avaDir,
+    statsMeasure="max",
+    flagShow=flagShow,
+    flagHue=True,
+    layer=cfgStats["GENERAL"]["layer"],
+)
+
+log.info("Plots have been saved to: %s" % outDir)
diff --git a/avaframe/tests/test_MoTUtils.py b/avaframe/tests/test_MoTUtils.py
index 373a0743c..f85a304b6 100644
--- a/avaframe/tests/test_MoTUtils.py
+++ b/avaframe/tests/test_MoTUtils.py
@@ -575,49 +575,6 @@ def test_setVariableForestParameters_noForest(tmp_path):
     assert result["File names"]["Tree diameter filename"] == "-"
 
 
-def test_MoTGenerateConfigs(tmp_path):
-    """Test MoTGenerateConfigs function"""
-    import configparser
-    import sys
-    from unittest.mock import MagicMock
-
-    # Create a mock module
-    mockModule = MagicMock()
-    mockModule.__name__ = "avaframe.com9MoTVoellmy.com9MoTVoellmy"
-
-    # Setup config
-    cfgMain = configparser.ConfigParser()
-    cfgMain["MAIN"] = {
-        "avalancheDir": str(tmp_path / "avaTest"),
-        "simTypeList": "null"
-    }
-
-    cfgInfo = None
-
-    # Mock the com1DFA functions
-    mockSimDict = {
-        "sim001": {
-            "cfgSim": configparser.ConfigParser()
-        }
-    }
-    mockInputSimFiles = {
-        "demFile": tmp_path / "dem.asc",
-        "relFiles": []
-    }
-
-    with patch('avaframe.in3Utils.MoTUtils.com1DFA.com1DFAPreprocess') as mockPreprocess:
-
-        mockPreprocess.return_value = (mockSimDict, tmp_path / "out", mockInputSimFiles, None)
-
-        # Call function
-        simDict, inputSimFiles = MoTUtils.MoTGenerateConfigs(cfgMain, cfgInfo, mockModule)
-
-        # Verify results
-        assert simDict == mockSimDict
-        assert inputSimFiles == mockInputSimFiles
-        mockPreprocess.assert_called_once()
-
-
 def test_RunAndCheckMoT_HighTimeStepCount():
     """Test that time step counter logs after 100 steps"""
     # Create output with 105 lines containing "Step" to trigger printCounter > 100
diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py
index aa2e0a944..469540d76 100644
--- a/avaframe/tests/test_com1DFA.py
+++ b/avaframe/tests/test_com1DFA.py
@@ -3554,17 +3554,40 @@ def test_com1DFAMainWithPathCfgInfo(tmp_path, caplog):
     # Pass pathlib.Path directly as cfgInfo
     cfgInfoPath = pathlib.Path(cfgDir)
 
-    with patch("avaframe.com1DFA.com1DFA.com1DFATools.createSimDictFromCfgs") as mockCreateSimDict:
+    with patch("avaframe.com1DFA.com1DFA.com1DFAPreprocess") as mockPreprocess:
 
-        mockCreateSimDict.return_value = ({}, {}, None, tmp_path / "out")
+        mockPreprocess.return_value = ({}, tmp_path / "out", {}, None)
 
-        # Call with directory Path - should route to createSimDictFromCfgs
+        # Call with directory Path - should route through com1DFAPreprocess
         try:
             com1DFA.com1DFAMain(cfgMain, cfgInfo=cfgInfoPath)
         except Exception:
             pass  # We expect it to fail due to missing simulations, but that's OK
 
-        # createSimDictFromCfgs should be called with the Path
-        mockCreateSimDict.assert_called_once()
-        callArgs = mockCreateSimDict.call_args
+        # com1DFAPreprocess should be called with the Path
+        mockPreprocess.assert_called_once()
+        callArgs = mockPreprocess.call_args
         assert callArgs[0][1] == cfgInfoPath
+
+
+def test_com1DFAPreprocessWithDirectoryPath(tmp_path):
+    """Test that com1DFAPreprocess handles a directory Path (batch mode)
+
+    When cfgInfo is a pathlib.Path pointing to a directory, com1DFAPreprocess should
+    route to createSimDictFromCfgs and return values in its standard order:
+    (simDict, outDir, inputSimFiles, simDFExisting)
+    """
+
+    dirPath = pathlib.Path(__file__).parents[0]
+    testPath = dirPath / "data" / "com1DFAConfigs"
+    inputDir = dirPath / "data" / "testCom1DFA2"
+    avaDir = pathlib.Path(tmp_path, "testCom1DFA")
+    shutil.copytree(inputDir, avaDir)
+    cfgMain = configparser.ConfigParser()
+    cfgMain["MAIN"] = {"avalancheDir": avaDir}
+
+    simDict, outDir, inputSimFiles, simDFExisting = com1DFA.com1DFAPreprocess(cfgMain, testPath)
+
+    assert len(simDict) == 16
+    assert simDFExisting is None
+    assert "demFile" in inputSimFiles
diff --git a/avaframe/tests/test_fileHandlerUtils.py b/avaframe/tests/test_fileHandlerUtils.py
index 6fbd75e27..210df49cf 100644
--- a/avaframe/tests/test_fileHandlerUtils.py
+++ b/avaframe/tests/test_fileHandlerUtils.py
@@ -494,3 +494,38 @@ def test_makeSimFromResDF_multiLayer_simName_reconstructed(tmp_path):
     assert "L1" not in dataDF["simName"].iloc[0]
     assert "L2" not in dataDF["simName"].iloc[0]
     assert dataDF["simName"].iloc[0] == simBase
+
+
+# --- Multi-layer makeSimDF tests ---
+
+
+def test_makeSimDF_layerColumn(tmp_path):
+    """makeSimDF should include a layer column populated from parseSimName"""
+    peakDir = tmp_path / "peakFiles"
+    peakDir.mkdir()
+
+    _createTestRaster(peakDir / "release1_abc123_null_psa_L1_ppr.asc")
+    _createTestRaster(peakDir / "release1_abc123_null_psa_L2_ppr.asc")
+    _createTestRaster(peakDir / "release1_abc123_null_psa_L1_pfv.asc")
+
+    dataDF = fU.makeSimDF(str(peakDir))
+
+    assert "layer" in dataDF.columns
+    # 3 files -> 3 rows
+    assert len(dataDF) == 3
+    layers = sorted(dataDF["layer"].tolist())
+    assert layers == ["L1", "L1", "L2"]
+
+
+def test_makeSimDF_layerColumn_singleLayer(tmp_path):
+    """makeSimDF layer column should be empty string for single-layer files"""
+    peakDir = tmp_path / "peakFiles"
+    peakDir.mkdir()
+
+    _createTestRaster(peakDir / "release1_abc123_null_dfa_ppr.asc")
+    _createTestRaster(peakDir / "release1_abc123_null_dfa_pfv.asc")
+
+    dataDF = fU.makeSimDF(str(peakDir))
+
+    assert "layer" in dataDF.columns
+    assert all(v == "" for v in dataDF["layer"].tolist())
diff --git a/avaframe/tests/test_getStats.py b/avaframe/tests/test_getStats.py
index cb61c082c..fbe3c2add 100644
--- a/avaframe/tests/test_getStats.py
+++ b/avaframe/tests/test_getStats.py
@@ -54,11 +54,20 @@ def test_getStats(tmp_path):
     # parameter dictionary
     varPar = 'relTh'
     inputDir = avaDirPeakFiles
-    peakValues = getStats.extractMaxValues(inputDir, avaDirtmp, varPar, restrictType='ppr', nameScenario='', parametersDict='')
+    peakValues, hasMultiLayer = getStats.extractMaxValues(
+        inputDir, avaDirtmp, varPar, restrictType="ppr", nameScenario="", parametersDict=""
+    )
 
     # call function to be tested test 2
     parametersDict = {'relTh': 1.0}
-    peakValues2 = getStats.extractMaxValues(inputDir, avaDirtmp, varPar, restrictType='', nameScenario='releaseScenario', parametersDict=parametersDict)
+    peakValues2, hasMultiLayer = getStats.extractMaxValues(
+        inputDir,
+        avaDirtmp,
+        varPar,
+        restrictType="",
+        nameScenario="releaseScenario",
+        parametersDict=parametersDict,
+    )
 
     assert peakValues['release1_null_dfa_1000000']['varPar'] == 1.0
     assert peakValues['release1_null_dfa_2000000']['varPar'] == 2.0
@@ -69,3 +78,66 @@ def test_getStats(tmp_path):
     assert peakValues['release1_null_dfa_2000000']['ppr']['min'] == 1.0
     assert peakValues['release1_null_dfa_2000000']['ppr']['mean'] == 3.0
     assert peakValues2['release1_null_dfa_1000000']['scenario'] == 'release1'
+    assert hasMultiLayer is False
+
+    with pytest.raises(ValueError) as e:
+        peakValues3, hasMultiLayer4 = getStats.extractMaxValues(
+            avaDirPeakFiles,
+            avaDirtmp,
+            varPar,
+            restrictType="",
+            nameScenario="",
+            parametersDict="",
+            layer="L1",
+        )
+    assert ("No multi-layer results detected but layer is specified") in str(e.value)
+
+
+def test_getStats_multiLayer(tmp_path):
+    """extractMaxValues should key by resType_layer for multi-layer files (e.g. ppr_l1, ppr_l2)"""
+
+    avaName = "avaStatsML"
+    avaDirtmp = pathlib.Path(tmp_path, avaName)
+    avaDirPeakFiles = avaDirtmp / "Outputs" / "com1DFA" / "peakFiles"
+    avaDirConfigFiles = avaDirtmp / "Outputs" / "com1DFA" / "configurationFiles"
+    os.makedirs(avaDirPeakFiles)
+    os.makedirs(avaDirConfigFiles)
+
+    # set path to existing test data files
+    dirPath = pathlib.Path(__file__).parents[0]
+    testDataDir = pathlib.Path(dirPath, "data")
+    data1 = testDataDir / "testGetStats_1000000_ppr.asc"
+    data2 = testDataDir / "testGetStats_2000000_ppr.asc"
+    cfg1 = testDataDir / "testGetStats_1000000.ini"
+
+    # create multi-layer peak files: L1 uses data1 values, L2 uses data2 values
+    shutil.copy(data1, avaDirPeakFiles / "release1_null_dfa_1000000_L1_ppr.asc")
+    shutil.copy(data2, avaDirPeakFiles / "release1_null_dfa_1000000_L2_ppr.asc")
+    shutil.copy(cfg1, avaDirConfigFiles / "release1_null_dfa_1000000.ini")
+
+    varPar = "relTh"
+    peakValues, hasMultiLayer = getStats.extractMaxValues(
+        avaDirPeakFiles, avaDirtmp, varPar, restrictType="", nameScenario="", parametersDict="", layer="L1"
+    )
+
+    simName = "release1_null_dfa_1000000"
+    # L1 and L2 should be separate keys, not overwriting each other
+    assert "ppr_l1" in peakValues[simName]
+    assert "ppr_l2" in peakValues[simName]
+    # L1 has data1 values (max=4), L2 has data2 values (max=10)
+    assert peakValues[simName]["ppr_l1"]["max"] == 4.0
+    assert peakValues[simName]["ppr_l2"]["max"] == 10.0
+    assert hasMultiLayer
+
+    with pytest.raises(ValueError) as e:
+        peakValues2, hasMultiLayer2 = getStats.extractMaxValues(
+            avaDirPeakFiles,
+            avaDirtmp,
+            varPar,
+            restrictType="",
+            nameScenario="",
+            parametersDict="",
+            layer="",
+        )
+    assert ("Multi-layer results detected but no layer specified.") in str(e.value)
+
diff --git a/avaframe/tests/test_probAna.py b/avaframe/tests/test_probAna.py
index 272cc16ae..a06506e4f 100644
--- a/avaframe/tests/test_probAna.py
+++ b/avaframe/tests/test_probAna.py
@@ -1239,3 +1239,118 @@ def test_createSample_invalid_method():
     # Check if appropriate error is raised
     with pytest.raises(AssertionError):
         createSample(testConfig, testParList)
+
+
+def _writeSmallAsc(filePath, nRows, nCols, data):
+    """Helper to write a minimal .asc raster file for testing"""
+    header = (
+        "ncols %d\n"
+        "nrows %d\n"
+        "xllcenter 0.0\n"
+        "yllcenter 0.0\n"
+        "cellsize 5.0\n"
+        "nodata_value -9999\n" % (nCols, nRows)
+    )
+    with open(filePath, "w") as f:
+        f.write(header)
+        for row in data:
+            f.write(" ".join(str(v) for v in row) + "\n")
+
+
+def test_probAnalysis_singleLayer_unchanged(tmp_path):
+    """Single-layer results with empty layer config work as before"""
+
+    # Create minimal directory structure
+    avaDir = tmp_path / "avaTest"
+    peakDir = avaDir / "Outputs" / "testMod" / "peakFiles"
+    peakDir.mkdir(parents=True)
+    outDir = avaDir / "Outputs" / "ana4Stats"
+    outDir.mkdir(parents=True)
+
+    # Create two single-layer ppr files with values > threshold(1.0) in different cells
+    nRows, nCols = 3, 3
+    data1 = [[0, 0, 0], [0, 5.0, 0], [0, 0, 0]]
+    data2 = [[0, 0, 0], [0, 5.0, 0], [0, 0, 2.0]]
+
+    _writeSmallAsc(peakDir / "rel1_abc123_com1_C_null_dfa_ppr.asc", nRows, nCols, data1)
+    _writeSmallAsc(peakDir / "rel1_def456_com1_C_null_dfa_ppr.asc", nRows, nCols, data2)
+
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {"peakLim": "1.0", "peakVar": "ppr", "layer": ""}
+
+    # Use a non-filtering modName so all files are included
+    analysisPerformed, contourDict = pA.probAnalysis(avaDir, cfg, "testMod")
+
+    assert analysisPerformed is True
+    # Both sims have center cell > threshold, so probability = 1.0 there
+    # Only sim2 has bottom-right > threshold, so probability = 0.5 there
+    probFile = outDir / "avaTest_prob__ppr_lim1.0.asc"
+    assert probFile.exists()
+    probData = np.loadtxt(probFile, skiprows=6)
+    assert probData.shape == (nRows, nCols)
+    # Center cell: both sims exceed threshold → 1.0
+    assert np.isclose(probData[1, 1], 1.0)
+    # Bottom-right: only 1 of 2 sims exceeds → 0.5
+    assert np.isclose(probData[2, 2], 0.5)
+
+
+def test_probAnalysis_multiLayer_withLayer(tmp_path):
+    """Multi-layer results with layer set filters to only the specified layer"""
+
+    avaDir = tmp_path / "avaTest"
+    peakDir = avaDir / "Outputs" / "com8MoTPSA" / "peakFiles"
+    peakDir.mkdir(parents=True)
+    outDir = avaDir / "Outputs" / "ana4Stats"
+    outDir.mkdir(parents=True)
+
+    nRows, nCols = 3, 3
+    # L1 files: high values everywhere
+    dataL1 = [[5.0, 5.0, 5.0], [5.0, 5.0, 5.0], [5.0, 5.0, 5.0]]
+    # L2 files: only center cell has high value
+    dataL2 = [[0, 0, 0], [0, 3.0, 0], [0, 0, 0]]
+
+    _writeSmallAsc(peakDir / "rel1_abc123_com8_C_null_psa_L1_ppr.asc", nRows, nCols, dataL1)
+    _writeSmallAsc(peakDir / "rel1_abc123_com8_C_null_psa_L2_ppr.asc", nRows, nCols, dataL2)
+    _writeSmallAsc(peakDir / "rel1_def456_com8_C_null_psa_L1_ppr.asc", nRows, nCols, dataL1)
+    _writeSmallAsc(peakDir / "rel1_def456_com8_C_null_psa_L2_ppr.asc", nRows, nCols, dataL2)
+
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {"peakLim": "1.0", "peakVar": "ppr", "layer": "L2"}
+
+    analysisPerformed, contourDict = pA.probAnalysis(
+        avaDir, cfg, "testMod", inputDir=avaDir / "Outputs" / "com8MoTPSA"
+    )
+
+    assert analysisPerformed is True
+    # Only L2 files should be included: both have center=3.0 > 1.0, rest=0
+    probFile = outDir / "avaTest_prob__ppr_L2_lim1.0.asc"
+    assert probFile.exists()
+    probData = np.loadtxt(probFile, skiprows=6)
+    # Center cell: both L2 sims exceed → 1.0
+    assert np.isclose(probData[1, 1], 1.0)
+    # Corner cells: no L2 sim exceeds → 0.0
+    assert np.isclose(probData[0, 0], 0.0)
+
+
+def test_probAnalysis_multiLayer_noLayer_errors(tmp_path):
+    """Multi-layer results without layer config raises an error"""
+
+    avaDir = tmp_path / "avaTest"
+    peakDir = avaDir / "Outputs" / "com8MoTPSA" / "peakFiles"
+    peakDir.mkdir(parents=True)
+    outDir = avaDir / "Outputs" / "ana4Stats"
+    outDir.mkdir(parents=True)
+
+    nRows, nCols = 3, 3
+    data = [[0, 0, 0], [0, 5.0, 0], [0, 0, 0]]
+
+    _writeSmallAsc(peakDir / "rel1_abc123_com8_C_null_psa_L1_ppr.asc", nRows, nCols, data)
+    _writeSmallAsc(peakDir / "rel1_abc123_com8_C_null_psa_L2_ppr.asc", nRows, nCols, data)
+
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {"peakLim": "1.0", "peakVar": "ppr", "layer": ""}
+
+    with pytest.raises(ValueError, match="Multi-layer results detected"):
+        pA.probAnalysis(
+            avaDir, cfg, "testMod", inputDir=avaDir / "Outputs" / "com8MoTPSA"
+        )
diff --git a/docs/_cfgFiles/ana3AIMECCfg.ini b/docs/_cfgFiles/ana3AIMECCfg.ini
index 2d6307177..aa94e890b 100644
--- a/docs/_cfgFiles/ana3AIMECCfg.ini
+++ b/docs/_cfgFiles/ana3AIMECCfg.ini
@@ -24,6 +24,10 @@ resTypes = ppr|pft|pfv
 # data result type for runout analysis (ppr, pft, pfv)
 runoutResType = ppr
 
+# layer to use for analysis (e.g. L1, L2). Leave empty for single-layer modules.
+# Required when analyzing multi-layer results (e.g. com8MoTPSA).
+runoutLayer =
+
 # limit value for evaluation of runout (depends on the runoutResType chosen)
 thresholdValue = 1
 
diff --git a/docs/moduleAna3AIMEC.rst b/docs/moduleAna3AIMEC.rst
index f90fce799..0befd4e90 100644
--- a/docs/moduleAna3AIMEC.rst
+++ b/docs/moduleAna3AIMEC.rst
@@ -107,6 +107,11 @@ To run
    - D - *modelType*: can be any descriptive string of the employed model (here dfa for dense flow  avalanche) 
    - E - *result type*: is pft (peak flow thickness) and pfv (peak flow velocity)
 
+   For multi-layer modules (e.g. com8MoTPSA), an additional layer component is inserted
+   before the result type: *A_B_C_D_L_E.*, where L is the layer identifier (e.g. L1, L2).
+   When analyzing multi-layer results, set ``runoutLayer`` in the AIMEC configuration
+   to select which layer to use for the runout analysis (e.g. ``runoutLayer = L1``).
+
 Theory
 -----------
 
diff --git a/docs/moduleAna4Stats.rst b/docs/moduleAna4Stats.rst
index 290583f1d..22d0b2bb4 100644
--- a/docs/moduleAna4Stats.rst
+++ b/docs/moduleAna4Stats.rst
@@ -77,6 +77,9 @@ also filtering is implemented. For other modules the only requirement is that si
 * *modelType*: for example **dfa** (dense flow avalanche) or **psa** (powder snow avalanche)
 * *resultType*: result variable, for example **pft**, **ppr**, **pfv**
 
+For multi-layer modules (e.g. com8MoTPSA), an additional layer component (e.g. **L1**, **L2**) is inserted
+before the resultType: ``..._modelType_layer_resultType.*``
+
 In order to run:
 
 * first go to ``AvaFrame/avaframe``