diff --git a/avaframe/ana1Tests/damBreak.py b/avaframe/ana1Tests/damBreakTest.py
similarity index 90%
rename from avaframe/ana1Tests/damBreak.py
rename to avaframe/ana1Tests/damBreakTest.py
index 15e0cab82..b5ff5f721 100644
--- a/avaframe/ana1Tests/damBreak.py
+++ b/avaframe/ana1Tests/damBreakTest.py
@@ -21,7 +21,7 @@
 log = logging.getLogger(__name__)
 
 
-def damBreakSol(avaDir, cfg, cfgC, outDirTest):
+def damBreakSol(avaDir, cfg, damBreakCfg, com1DFACfg, outDirTest):
     """ Compute analytical Flow thickness and velocity for dam break test case
 
         for a granular flow over a dry rough sloping bed with the Savage Hutter model
@@ -32,8 +32,10 @@ def damBreakSol(avaDir, cfg, cfgC, outDirTest):
             path to avalanche directory
         cfg: configParser object
             main avaframe configuration - here showPlot flag used
-        cfgC: configParser object
-            configuration setting for avalanche simulation including DAMBREAK section
+        damBreakCfg: configParser object
+            dam break configuration
+        com1DFACfg: configParser object
+            com1DFA configuration
         outDirTest: pathlib path
             path to output directory
 
@@ -58,21 +60,22 @@ def damBreakSol(avaDir, cfg, cfgC, outDirTest):
 
     # Set Parameters
     # Coordinate system chosen in the direction of the inclined plane
-    g = cfgC['GENERAL'].getfloat('gravAcc')       # acceleration due to gravity [ms-2]
-    phiDeg = cfgC['DAMBREAK'].getfloat('phi')       # slope angle [°]
-    deltaDeg = cfgC['DAMBREAK'].getfloat('delta')       # friction angle [°]
+    cfgGen = com1DFACfg['GENERAL']
+    g = cfgGen.getfloat('gravAcc')       # acceleration due to gravity [ms-2]
+    phiDeg = damBreakCfg.getfloat('phi')       # slope angle [°]
+    deltaDeg = damBreakCfg.getfloat('delta')       # friction angle [°]
     phi = np.radians(phiDeg)                          # slope angle [rad]
     delta = np.radians(deltaDeg)                        # bed friction angle [rad]
     gz = g * np.cos(phi)                          # projection of g perpendicular to the inclined plane
     m0 = gz * (np.tan(phi) - np.tan(delta))       # constant x-acceleration resulting from gravity and friction force
-    hL = cfgC['GENERAL'].getfloat('relTh')        # initial height [m] in Riemann problem in state 1 (x<0), hR (x>0)=0
+    hL = cfgGen.getfloat('relTh')        # initial height [m] in Riemann problem in state 1 (x<0), hR (x>0)=0
     cL = np.sqrt(gz * hL)
-    x0R = cfgC['DAMBREAK'].getfloat('xBack') / np.cos(phi)
-    tSave = cfgC['DAMBREAK'].getfloat('tSave')
+    x0R = damBreakCfg.getfloat('xBack') / np.cos(phi)
+    tSave = damBreakCfg.getfloat('tSave')
     # Define time [0-1] seconds and space [-2,2] meters domains multiplied times 100
-    tAna = np.arange(0, cfgC['DAMBREAK'].getfloat('tEnd'), cfgC['DAMBREAK'].getfloat('dt'))
-    x = np.arange(cfgC['DAMBREAK'].getfloat('xStart') / np.cos(phi) , cfgC['DAMBREAK'].getfloat('xEnd') / np.cos(phi),
-                  cfgC['DAMBREAK'].getfloat('dx'))
+    tAna = np.arange(0, damBreakCfg.getfloat('tEnd'), damBreakCfg.getfloat('dt'))
+    x = np.arange(damBreakCfg.getfloat('xStart') / np.cos(phi), damBreakCfg.getfloat('xEnd') / np.cos(phi),
+                  damBreakCfg.getfloat('dx'))
     nt = len(tAna)
     nx = len(x)
     # Initialise Flow thickness solution and velocity
@@ -227,13 +230,19 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solDam, fieldHeader, cfg,
     xAna = solDam['xAna']
     hAna = solDam['hAna']
     uAna = solDam['uAna']
-    # get plots limits
-    limits = outAna1Plots.getPlotLimits(cfgDam, fieldsList, fieldHeader)
     hErrorL2Array = np.zeros((len(fieldsList)))
     vhErrorL2Array = np.zeros((len(fieldsList)))
     hErrorLMaxArray = np.zeros((len(fieldsList)))
     vhErrorLMaxArray = np.zeros((len(fieldsList)))
     count = 0
+    tSave = cfgDam.getfloat('tSave')
+    indT = min(np.searchsorted(timeList, tSave), min(len(timeList)-1, len(fieldsList)-1))
+    tSave = timeList[indT]
+    # get plots limits
+    if cfgDam.getboolean('plotSequence'):
+        limits = outAna1Plots.getPlotLimits(cfgDam, fieldsList, fieldHeader)
+    else:
+        limits = outAna1Plots.getPlotLimits(cfgDam, [fieldsList[indT]], fieldHeader)
     # run the comparison routine for each saved time step
     for t, field in zip(timeList, fieldsList):
         # get similartiy solution h, u at required time step
@@ -251,9 +260,9 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solDam, fieldHeader, cfg,
 
         # make analytical results 2D
         nrows, ncols = np.shape(hNumPlus)
-        O = np.ones((nrows, ncols))
-        hDamPlus = O*hDamPlus
-        uDamPlus = O*uDamPlus
+        OnesRaster = np.ones((nrows, ncols))
+        hDamPlus = OnesRaster*hDamPlus
+        uDamPlus = OnesRaster*uDamPlus
 
         hEL2Plus, hL2RelPlus, hLmaxPlus, hLmaxRelPlus = anaTools.normL2Scal(hDamPlus, hNumPlus, cellSize,
                                                                             np.cos(phiRad))
@@ -273,15 +282,12 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solDam, fieldHeader, cfg,
         log.debug("L2 %s error on the Flow thickness at t=%.2f s is : %.4f" % (title, t, hEL2Plus))
         log.debug("L2 %s error on the momentum at t=%.2f s is : %.4f" % (title, t, vhL2Plus))
         # Make all individual time step comparison plot
-        if cfgDam.getboolean('plotSequence'):
+        if cfgDam.getboolean('plotSequence') or t == tSave:
             outAna1Plots.plotComparisonDam(cfgDam, simHash, fieldsList[0], field, fieldHeader, solDam, t, limits,
                                            outDirTest)
         count = count + 1
 
-    tSave = cfgDam.getfloat('tSave')
-    indT = min(np.searchsorted(timeList, tSave), min(len(timeList)-1, len(fieldsList)-1))
-    tSave = timeList[indT]
-    if cfgDam.getboolean('plotErrorTime') and len(timeList)>1:
+    if cfgDam.getboolean('plotErrorTime') and len(timeList) > 1:
         # Create result plots
 
         outAna1Plots.plotErrorTime(timeList, hErrorL2Array, hErrorLMaxArray, vhErrorL2Array, vhErrorLMaxArray,
diff --git a/avaframe/ana1Tests/damBreakTestCfg.ini b/avaframe/ana1Tests/damBreakTestCfg.ini
new file mode 100644
index 000000000..c9e1f3e07
--- /dev/null
+++ b/avaframe/ana1Tests/damBreakTestCfg.ini
@@ -0,0 +1,140 @@
+### Config File - This file contains the main settings for the dam break test module
+## Set your parameters
+# This file will be overridden by local_damBreakTestCfg.ini if it exists
+# So copy this file to local_damBreakTestCfg.ini, adjust your variables there
+
+[DAMBREAK]
+#sphKernel radius for the com1DFA simulations
+sphKernelRadius = 6|5|4|3
+
+# slope angle [°]
+phi = 22
+# bed friction angle [°]
+delta = 21
+u0 = 0.
+# time end and resolution for the analytic solution
+tEnd = 30
+dt = 0.1
+# space resolution
+dx = 0.5
+# domain extend for error computation
+# in x direction
+xStart = -200
+xEnd = 220
+# in y direction
+yStart = -50
+yEnd = 50
+# start x position of the dam
+xBack = -120
+# xFront = 0
+
+# set time step of analysis
+tSave = 15
+# use only the component of the velocity/momentum in the flow direction (vector going down the inclined plane in x dir)
+projectVelocity = False
+
+#++++Plotting parameters++++++
+# list of parameters to display in the summary plot (list of parameters separated by |)
+paramInfo = sphKernelRadius|aPPK|nPPK0|nPart
+# list of parameters to display in the convergence plot (list of parameters separated by |)
+paramConvInfo = sphKR0|nPPK0|cMax
+# plotting flags
+# only analyze and plot the tSave time step
+onlyLast = False
+# plot error function of time for each simulation
+plotErrorTime = False
+# plot individual figure for the h, hu and u error for each saved time step
+plotSequence = False
+# use relative difference
+relativError = True
+# when plotting, the domain extent is scaleCoef times bigger than the data extent
+scaleCoef = 1.05
+
+
+
+[com1DFA_override]
+# use default com1DFA config as base configuration (True) and override following parameters
+# if False and local is available use local
+defaultConfig = True
+#++++++++++++++++ Simulation type
+# list of simulations that shall be performed (null, ent, res, entres, available (use all available input data))
+simTypeList = null
+#+++++++++++++ Output++++++++++++
+# desired result Parameters (ppr, pfd, pfv, FD, FV, P, particles) - separated by |
+resType = FT|FV|Vx|Vy|Vz
+# saving time step, i.e.  time in seconds (first and last time step are always saved)
+# 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)
+# option 2: explicitly list all desired time steps (closest to actual computational time step) separated by | (example tSteps = 1|50.2|100)
+# NOTE: initial and last time step are always saved!
+tSteps = 0:1
+
+#++++++++++++++++ particle Initialisation +++++++++
+# particle initialisation step - set iniStep to True to activate redistribution of particles to reduce SPH force
+# this is in a development stage - hence parameters are set for development and will be adjusted after extensive testing
+iniStep = True
+# max number of iterations - high number might cause significant increase in computational time
+maxIterations = 30
+# buffer zone factor multiplied with sphKernelRadius
+bufferZoneFactor = 4
+
+#+++++++++SNOW properties
+# True if release thickness should be read from shapefile file; if False - relTh read from ini file
+relThFromShp = False
+relTh = 1
+
+#++++++++++++Time stepping parameters
+# End time [s]
+tEnd = 20
+# to use a variable time step (time step depends on kernel radius)
+sphKernelRadiusTimeStepping = True
+# courant number
+cMax = 0.01
+
+#+++++++++++++SPH parameters
+# SPH gradient option
+# 0) No pressure gradients computed
+# 1) SamosAT style (no reprojecion on the surface, dz = 0 and gz is used)
+# 2) SamostAt but done in the cartesian coord system (will hopefully allow us to add the earth pressure coef)
+# 3) SamostAt but done in the local coord system (will hopefully allow us to add the earth pressure coef)
+# and this time reprojecion on the surface, dz not 0 and g3 is used
+sphOption = 2
+# sph kernel smoothing length [m]
+sphKernelRadius = 3
+# Choice of artificial viscosity
+# 0) No artificial viscosity
+# 1) SAMOS artificial viscosity
+# 2) Ata artificial viscosity
+viscOption = 0
+
+#++++++++++++++++ Particles
+# number of particles defined by: MPPDIR= mass per particle direct, MPPDH= mass per particles through release thickness,
+# MPPKR= mass per particles through number of particles per kernel radius
+massPerParticleDeterminationMethod = MPPKR
+# number of particles per kernel radius (if MPPKR is used)
+# is computed with: nPPK = nPPK0 * (sphKR/sphKR0)^aPPK
+# where sphKR is the sphKernelRadius specified further up
+# reference kernel radius
+sphKR0 = 5
+# reference number of particles per kernel radius
+nPPK0 = 15
+# variation of nppK exponent
+aPPK = 0|-0.5|-1|-2|-3
+
+# if splitOption=0
+# threshold for splitting particles, split if: mPart > massPerPart x thresholdMassSplit
+thresholdMassSplit = 5
+
+
+#+++++++++++++Mesh and interpolation
+# remesh the input DEM
+# expected mesh size [m]
+meshCellSize = 3
+
+#++++++++++++Friction model
+# add the friction using an explicit formulation (1)
+# 0 use an implicit method
+explicitFriction = 1
+# friction type (samosAT, Coulomb)
+frictModel = Coulomb
+#+++++++++++++SamosAt friction model
+mu = 0.3838640350354158
diff --git a/avaframe/ana1Tests/energyLineTest.py b/avaframe/ana1Tests/energyLineTest.py
index 72e103ee4..8b1035700 100644
--- a/avaframe/ana1Tests/energyLineTest.py
+++ b/avaframe/ana1Tests/energyLineTest.py
@@ -135,11 +135,12 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
     zMin = zLim[0]
     ax2.vlines(x=avaProfileMass['s'][-1], ymin=zMin, ymax=avaProfileMass['z'][-1],
                color='r', linestyle='--')
-    ax2.vlines(x=sIntersection, color='b', ymin=zMin, ymax=zIntersection, linestyle='--')
+    if ~np.isnan(sIntersection):
+        ax2.vlines(x=sIntersection, color='b', ymin=zMin, ymax=zIntersection, linestyle='--')
+        ax2.hlines(y=zIntersection, color='b', xmin=0, xmax=sIntersection, linestyle='--')
 
     ax2.hlines(y=avaProfileMass['z'][-1], xmin=0, xmax=avaProfileMass['s'][-1],
                color='r', linestyle='--')
-    ax2.hlines(y=zIntersection, color='b', xmin=0, xmax=sIntersection, linestyle='--')
 
     ax2.set_xlabel('s [m]')
     ax2.set_ylabel('z [m]')
@@ -183,8 +184,8 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
     rmseVelocityElevation = resultEnergyTest['rmseVelocityElevation']
     errorS = abs(runOutSError)
     errorZ = abs(runOutZError)
-    sMin = min(avaProfileMass['s'][-1], sIntersection) - max(errorS, 0)
-    sMax = max(avaProfileMass['s'][-1], sIntersection) + max(errorS, 0)
+    sMin = np.nanmin(np.array([avaProfileMass['s'][-1], sIntersection])) - np.nanmax(np.array([errorS, 0]))
+    sMax = np.nanmax(np.array([avaProfileMass['s'][-1], sIntersection])) + np.nanmax(np.array([errorS, 0]))
     zMin = avaProfileMass['z'][-1] + min(slopeExt*(sMax - avaProfileMass['s'][-1]), 0 - 2*errorZ)-z0
     zMax = avaProfileMass['z'][0] - sMin*np.tan(min(runOutAngleRad, alphaRad))-z0
     if avaProfileMass['z'][-1] == zIntersection:
@@ -193,11 +194,12 @@ def generateCom1DFAEnergyPlot(avalancheDir, energyLineTestCfg, com1DFACfg, avaPr
         zMin = zMin - (zMax-zMin)*0.1
     ax3.vlines(x=avaProfileMass['s'][-1], ymin=zMin, ymax=avaProfileMass['z'][-1]-z0,
                color='r', linestyle='--')
-    ax3.vlines(x=sIntersection, color='b', ymin=zMin, ymax=zIntersection-z0, linestyle='--')
+    if ~np.isnan(sIntersection):
+        ax3.vlines(x=sIntersection, color='b', ymin=zMin, ymax=zIntersection-z0, linestyle='--')
+        ax3.hlines(y=zIntersection-z0, color='b', xmin=sMin, xmax=sIntersection, linestyle='--')
 
     ax3.hlines(y=avaProfileMass['z'][-1]-z0, xmin=sMin, xmax=avaProfileMass['s'][-1],
                color='r', linestyle='--')
-    ax3.hlines(y=zIntersection-z0, color='b', xmin=sMin, xmax=sIntersection, linestyle='--')
 
     ax3.set_xlabel('s [m]')
     ax3.set_ylabel('$\Delta z$ [m]')
@@ -310,8 +312,14 @@ def getAlphaProfileIntersection(energyLineTestCfg, avaProfileMass, mu, csz):
     zP1 = z[idx+1]
     zA0 = alphaLine[idx]
     zA1 = alphaLine[idx+1]
-    sIntersection = s0 + (s1-s0)*(zA0-zP0)/((zP1-zP0)-(zA1-zA0))
-    zIntersection = zP0 + (sIntersection-s0) * (zP1-zP0) / (s1-s0)
+    if s0 == 0:
+        sIntersection = np.nan
+        zIntersection = np.nan
+        s[-1] = avaProfileMass['s'][-1]
+        z[-1] = avaProfileMass['z'][-1]
+    else:
+        sIntersection = s0 + (s1-s0)*(zA0-zP0)/((zP1-zP0)-(zA1-zA0))
+        zIntersection = zP0 + (sIntersection-s0) * (zP1-zP0) / (s1-s0)
     if coefExt > 0:
         s[-1] = sIntersection
         z[-1] = zIntersection
@@ -363,7 +371,10 @@ def getEnergyInfo(avaProfileMass, g, mu, sIntersection, zIntersection, runOutAng
     GK = sIntersection * mu
     z_end = avaProfileMass['z'][0] - GK
     zEneTarget = avaProfileMass['z'][0] - avaProfileMass['s'] * mu
-    sGeomL = [avaProfileMass['s'][0], sIntersection]
+    if ~np.isnan(sIntersection):
+        sGeomL = [avaProfileMass['s'][0], sIntersection]
+    else:
+        sGeomL = [avaProfileMass['s'][0], avaProfileMass['s'][-1]]
     zGeomL = [avaProfileMass['z'][0], z_end]
     # compute errors
     # rmse of the energy height
diff --git a/avaframe/ana1Tests/energyLineTestCfg.ini b/avaframe/ana1Tests/energyLineTestCfg.ini
index a94c1a1a5..09a4b3db8 100644
--- a/avaframe/ana1Tests/energyLineTestCfg.ini
+++ b/avaframe/ana1Tests/energyLineTestCfg.ini
@@ -93,3 +93,5 @@ frictModel = Coulomb
 #+++++++++++++General Friction parameters
 # tan of bed friction angle used for: samosAT, Coulomb, Voellmy
 mu = 0.4
+#+++++++++++++Voellmy friction model
+xsi = 10000.
diff --git a/avaframe/ana1Tests/figAnalytic_energyLineTestCfg.ini b/avaframe/ana1Tests/figAnalytic_energyLineTestCfg.ini
new file mode 100644
index 000000000..66c68f646
--- /dev/null
+++ b/avaframe/ana1Tests/figAnalytic_energyLineTestCfg.ini
@@ -0,0 +1,82 @@
+### Config File - This file contains the main settings for the energy line test module
+## Set your parameters
+# This file will be overridden by local_energyLineTestCfg.ini if it exists
+# So copy this file to local_energyLineTestCfg.ini, adjust your variables there
+
+[energyLineTest]
+# do you want to run the DFA module (all results in the Outputs/com1DFA forlder will be deleted)
+runDFAModule = True
+# for extending the path at the bottom, use the last points of the mass average path to get the slope
+# (all points at a distance nCellsMinExtend x cellSize < distance from the end)
+nCellsExtrapolation = 4
+# plot the mass averaged points function of the corrected s too (only available if pathFromPart = True)
+plotScor = False
+
+# True to get the path from particles, from fields otherwise
+# if True, particles need to be saved in resType
+# if False, FT|FV|FM need to be saved in resType
+pathFromPart = True
+
+# for plotting the results, add the shifted alpha line
+shiftedAlphaLine = False
+
+[com1DFA_override]
+# use default com1DFA config as base configuration (True) and override following parameters
+# if False and local is available use local
+defaultConfig = True
+#++++++++++++++++ Simulation type
+# list of simulations that shall be performed (null, ent, res, entres, available (use all available input data))
+simTypeList = null
+
+#+++++++++++++ Output++++++++++++
+# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, Vx, Vy, Vz, TA, particles) - separated by |
+resType = FT|FV|FM|particles
+# saving time step, i.e.  time in seconds (first and last time step are always saved)
+# 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)
+# option 2: explicitly list all desired time steps (closest to actual computational time step) separated by | (example tSteps = 1|50.2|100)
+# NOTE: initial and last time step are always saved!
+tSteps = 0:1
+
+#+++++++++SNOW properties
+#+++++Release thickness++++
+# True if release thickness should be read from shapefile file; if False - relTh read from ini file
+relThFromShp = False
+# release thickness (only considered if relThFromShp=False)
+relTh = 1
+
+#++++++++++++Time stepping parameters
+# to use a variable time step (time step depends on kernel radius)
+sphKernelRadiusTimeStepping = True
+# Upper time step limit coefficient if option sphKernelRadiusTimeStepping is chosen.
+cMax = 0.02
+# stopCriterion (stops when massFlowing<0.01*peakMassFlowing or ke<0.01*pke)
+stopCrit = 0.0
+
+#+++++++++++++SPH parameters
+# SPH gradient option
+# it is possible to set other values here (2: SamosAT done in the cartesian coord system (reprojecion on the surface, dz != 0 and g3 is used))
+sphOption = 0
+# Choice of artificial viscosity
+# it is possible to set other values here  (1: add some artificial viscosity)
+viscOption = 0
+
+#++++++++++++++++ Particles
+# number of particles defined by: MPPDIR= mass per particle direct, MPPDH= mass per particles through release thickness,
+# MPPKR= mass per particles through number of particles per kernel radius
+massPerParticleDeterminationMethod = MPPKR
+
+#+++++++++++++Flow model parameters+++++++++
+# take curvature term into account in the gravity acceleration term for computing the friction force
+# (and gradient if curvAccInGradient=1)
+# 0 if deactivated, 1 if activated
+curvAccInFriction = 0
+
+#++++++++++++Friction model
+# add the friction using an explicit formulation (1)
+# 0 use an implicit method
+explicitFriction = 1
+# friction type (samosAT, Coulomb, Voellmy)
+frictModel = Coulomb
+#+++++++++++++General Friction parameters
+# tan of bed friction angle used for: samosAT, Coulomb, Voellmy
+mu = 0.4
diff --git a/avaframe/ana1Tests/figConvergence_damBreakTestCfg.ini b/avaframe/ana1Tests/figConvergence_damBreakTestCfg.ini
new file mode 100644
index 000000000..eb484b601
--- /dev/null
+++ b/avaframe/ana1Tests/figConvergence_damBreakTestCfg.ini
@@ -0,0 +1,50 @@
+### Config File - This file contains the main settings for the dam break test module
+## Set your parameters
+# This file will be overridden by local_damBreakTestCfg.ini if it exists
+# So copy this file to local_damBreakTestCfg.ini, adjust your variables there
+
+[DAMBREAK]
+#sphKernel radius for the com1DFA simulations
+sphKernelRadius = 6|5|4|3
+
+# set time step of analysis
+tSave = 15
+
+#++++Plotting parameters++++++
+# list of parameters to display in the summary plot (list of parameters separated by |)
+paramInfo = sphKernelRadius|aPPK|nPPK0|nPart
+# list of parameters to display in the convergence plot (list of parameters separated by |)
+paramConvInfo = sphKR0|nPPK0|cMax
+# plotting flags
+# only analyze and plot the tSave time step
+onlyLast = False
+# plot error function of time for each simulation
+plotErrorTime = False
+# plot individual figure for the h, hu and u error for each saved time step
+plotSequence = False
+# use relative difference
+relativError = True
+# when plotting, the domain extent is scaleCoef times bigger than the data extent
+scaleCoef = 1.05
+
+
+
+[com1DFA_override]
+# use default com1DFA config as base configuration (True) and override following parameters
+# if False and local is available use local
+defaultConfig = True
+
+#++++++++++++Time stepping parameters
+# End time [s]
+tEnd = 20
+# courant number
+# it is possible to set other values here (0.04|0.02|0.01|0.005)
+cMax = 0.02
+
+#++++++++++++++++ Particles
+# reference number of particles per kernel radius
+# it is possible to set other values here (20|30|40)
+nPPK0 = 15
+# variation of nppK exponent
+# it is possible to set other values here (0|-0.5|-1|-1.5|-2|-2.5|-3)
+aPPK = 0|-0.5|-1|-2|-3
diff --git a/avaframe/ana1Tests/figConvergence_simiSolTestCfg.ini b/avaframe/ana1Tests/figConvergence_simiSolTestCfg.ini
new file mode 100644
index 000000000..02144ae44
--- /dev/null
+++ b/avaframe/ana1Tests/figConvergence_simiSolTestCfg.ini
@@ -0,0 +1,50 @@
+### Config File - This file contains the main settings for the similarity solution test module
+## Set your parameters
+# This file will be overridden by local_simiSolTestCfg.ini if it exists
+# So copy this file to local_simiSolTestCfg.ini, adjust your variables there
+
+[SIMISOL]
+#sphKernel radius for the com1DFA simulations
+sphKernelRadius = 10|8|6|5|4|3
+
+# save analysis plots at time step dtSol
+dtSol = 50.
+
+# which time step should be saved
+tSave = 20
+
+#++++Plotting parameters++++++
+# list of parameters to display in the summary plot (list of parameters separated by |)
+paramInfo = sphKernelRadius|aPPK|nPPK0|nPart
+# list of parameters to display in the convergence plot (list of parameters separated by |)
+paramConvInfo = sphKR0|nPPK0|cMax
+# plotting flags
+# only analyze and plot the tSave time step
+onlyLast = False
+# plot error function of time for each simulation
+plotErrorTime = False
+# plot individual figure for the h, hu and u error for each saved time step
+plotSequence = False
+# use relative difference
+relativError = True
+# when plotting, the domain extent is scaleCoef times bigger than the data extent
+scaleCoef = 1.02
+
+
+[com1DFA_override]
+# use default com1DFA config as base configuration (True) and override following parameters
+# if False and local is available use local
+defaultConfig = True
+
+#++++++++++++Time stepping parameters
+# End time [s]
+tEnd = 30
+# it is possible to set other values here (0.04|0.02|0.005)
+cMax = 0.01
+
+# reference number of particles per kernel radius
+# it is possible to set other values here (20|30|40)
+nPPK0 = 15
+# variation of nppK exponent
+# it is possible to set other values here (0|-0.5|-1|-1.5|-2|-2.5|-3)
+aPPK = 0|-0.5|-1|-2|-3
diff --git a/avaframe/ana1Tests/figCurvature_energyLineTestCfg.ini b/avaframe/ana1Tests/figCurvature_energyLineTestCfg.ini
new file mode 100644
index 000000000..a726d5b5d
--- /dev/null
+++ b/avaframe/ana1Tests/figCurvature_energyLineTestCfg.ini
@@ -0,0 +1,82 @@
+### Config File - This file contains the main settings for the energy line test module
+## Set your parameters
+# This file will be overridden by local_energyLineTestCfg.ini if it exists
+# So copy this file to local_energyLineTestCfg.ini, adjust your variables there
+
+[energyLineTest]
+# do you want to run the DFA module (all results in the Outputs/com1DFA forlder will be deleted)
+runDFAModule = True
+# for extending the path at the bottom, use the last points of the mass average path to get the slope
+# (all points at a distance nCellsMinExtend x cellSize < distance from the end)
+nCellsExtrapolation = 4
+# plot the mass averaged points function of the corrected s too (only available if pathFromPart = True)
+plotScor = False
+
+# True to get the path from particles, from fields otherwise
+# if True, particles need to be saved in resType
+# if False, FT|FV|FM need to be saved in resType
+pathFromPart = True
+
+# for plotting the results, add the shifted alpha line
+shiftedAlphaLine = True
+
+[com1DFA_override]
+# use default com1DFA config as base configuration (True) and override following parameters
+# if False and local is available use local
+defaultConfig = True
+#++++++++++++++++ Simulation type
+# list of simulations that shall be performed (null, ent, res, entres, available (use all available input data))
+simTypeList = null
+
+#+++++++++++++ Output++++++++++++
+# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, Vx, Vy, Vz, TA, particles) - separated by |
+resType = FT|FV|FM|particles
+# saving time step, i.e.  time in seconds (first and last time step are always saved)
+# 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)
+# option 2: explicitly list all desired time steps (closest to actual computational time step) separated by | (example tSteps = 1|50.2|100)
+# NOTE: initial and last time step are always saved!
+tSteps = 0:1
+
+#+++++++++SNOW properties
+#+++++Release thickness++++
+# True if release thickness should be read from shapefile file; if False - relTh read from ini file
+relThFromShp = False
+# release thickness (only considered if relThFromShp=False)
+relTh = 1
+
+#++++++++++++Time stepping parameters
+# to use a variable time step (time step depends on kernel radius)
+sphKernelRadiusTimeStepping = True
+# Upper time step limit coefficient if option sphKernelRadiusTimeStepping is chosen.
+cMax = 0.02
+# stopCriterion (stops when massFlowing<0.01*peakMassFlowing or ke<0.01*pke)
+stopCrit = 0.0
+
+#+++++++++++++SPH parameters
+# SPH gradient option
+# it is possible to set other values here (2: SamosAT done in the cartesian coord system (reprojecion on the surface, dz != 0 and g3 is used))
+sphOption = 0
+# Choice of artificial viscosity
+# it is possible to set other values here  (1: add some artificial viscosity)
+viscOption = 0
+
+#++++++++++++++++ Particles
+# number of particles defined by: MPPDIR= mass per particle direct, MPPDH= mass per particles through release thickness,
+# MPPKR= mass per particles through number of particles per kernel radius
+massPerParticleDeterminationMethod = MPPKR
+
+#+++++++++++++Flow model parameters+++++++++
+# take curvature term into account in the gravity acceleration term for computing the friction force
+# (and gradient if curvAccInGradient=1)
+# 0 if deactivated, 1 if activated
+curvAccInFriction = 1
+
+#++++++++++++Friction model
+# add the friction using an explicit formulation (1)
+# 0 use an implicit method
+explicitFriction = 1
+# friction type (samosAT, Coulomb, Voellmy)
+frictModel = Coulomb
+#+++++++++++++General Friction parameters
+# tan of bed friction angle used for: samosAT, Coulomb, Voellmy
+mu = 0.4
diff --git a/avaframe/ana1Tests/figPressure_energyLineTestCfg.ini b/avaframe/ana1Tests/figPressure_energyLineTestCfg.ini
new file mode 100644
index 000000000..f7f5d27e6
--- /dev/null
+++ b/avaframe/ana1Tests/figPressure_energyLineTestCfg.ini
@@ -0,0 +1,82 @@
+### Config File - This file contains the main settings for the energy line test module
+## Set your parameters
+# This file will be overridden by local_energyLineTestCfg.ini if it exists
+# So copy this file to local_energyLineTestCfg.ini, adjust your variables there
+
+[energyLineTest]
+# do you want to run the DFA module (all results in the Outputs/com1DFA forlder will be deleted)
+runDFAModule = True
+# for extending the path at the bottom, use the last points of the mass average path to get the slope
+# (all points at a distance nCellsMinExtend x cellSize < distance from the end)
+nCellsExtrapolation = 4
+# plot the mass averaged points function of the corrected s too (only available if pathFromPart = True)
+plotScor = False
+
+# True to get the path from particles, from fields otherwise
+# if True, particles need to be saved in resType
+# if False, FT|FV|FM need to be saved in resType
+pathFromPart = True
+
+# for plotting the results, add the shifted alpha line
+shiftedAlphaLine = False
+
+[com1DFA_override]
+# use default com1DFA config as base configuration (True) and override following parameters
+# if False and local is available use local
+defaultConfig = True
+#++++++++++++++++ Simulation type
+# list of simulations that shall be performed (null, ent, res, entres, available (use all available input data))
+simTypeList = null
+
+#+++++++++++++ Output++++++++++++
+# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, Vx, Vy, Vz, TA, particles) - separated by |
+resType = FT|FV|FM|particles
+# saving time step, i.e.  time in seconds (first and last time step are always saved)
+# 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)
+# option 2: explicitly list all desired time steps (closest to actual computational time step) separated by | (example tSteps = 1|50.2|100)
+# NOTE: initial and last time step are always saved!
+tSteps = 0:1
+
+#+++++++++SNOW properties
+#+++++Release thickness++++
+# True if release thickness should be read from shapefile file; if False - relTh read from ini file
+relThFromShp = False
+# release thickness (only considered if relThFromShp=False)
+relTh = 1
+
+#++++++++++++Time stepping parameters
+# to use a variable time step (time step depends on kernel radius)
+sphKernelRadiusTimeStepping = True
+# Upper time step limit coefficient if option sphKernelRadiusTimeStepping is chosen.
+cMax = 0.02
+# stopCriterion (stops when massFlowing<0.01*peakMassFlowing or ke<0.01*pke)
+stopCrit = 0.0
+
+#+++++++++++++SPH parameters
+# SPH gradient option
+# it is possible to set other values here (2: SamosAT done in the cartesian coord system (reprojecion on the surface, dz != 0 and g3 is used))
+sphOption = 2
+# Choice of artificial viscosity
+# it is possible to set other values here  (1: add some artificial viscosity)
+viscOption = 1
+
+#++++++++++++++++ Particles
+# number of particles defined by: MPPDIR= mass per particle direct, MPPDH= mass per particles through release thickness,
+# MPPKR= mass per particles through number of particles per kernel radius
+massPerParticleDeterminationMethod = MPPKR
+
+#+++++++++++++Flow model parameters+++++++++
+# take curvature term into account in the gravity acceleration term for computing the friction force
+# (and gradient if curvAccInGradient=1)
+# 0 if deactivated, 1 if activated
+curvAccInFriction = 0
+
+#++++++++++++Friction model
+# add the friction using an explicit formulation (1)
+# 0 use an implicit method
+explicitFriction = 1
+# friction type (samosAT, Coulomb, Voellmy)
+frictModel = Coulomb
+#+++++++++++++General Friction parameters
+# tan of bed friction angle used for: samosAT, Coulomb, Voellmy
+mu = 0.4
diff --git a/avaframe/ana1Tests/rotationTestCfg.ini b/avaframe/ana1Tests/rotationTestCfg.ini
index b022c96d2..2f840df2d 100644
--- a/avaframe/ana1Tests/rotationTestCfg.ini
+++ b/avaframe/ana1Tests/rotationTestCfg.ini
@@ -7,7 +7,7 @@
 # for the report (which computational module do you use)
 comModule = com1DFA
 # do you want to run the DFA module (all results in the Outputs/com1DFA folder will be deleted)
-runDFAModule = False
+runDFAModule = True
 
 [energyLineTest_override]
 # use default energyLineTest config as base configuration (True) and override following parameters
@@ -51,6 +51,7 @@ relTh = 1
 #++++++++++++Time stepping parameters
 # to use a variable time step (time step depends on kernel radius)
 sphKernelRadiusTimeStepping = True
+# Upper time step limit coefficient if option sphKernelRadiusTimeStepping is chosen.
 cMax = 0.02
 # stopCriterion (stops when massFlowing<0.01*peakMassFlowing or ke<0.01*pke)
 stopCrit = 0.0
@@ -77,10 +78,15 @@ massPerParticleDeterminationMethod = MPPKR
 #+++++++++++++Flow model parameters+++++++++
 # subgridMixingFactor
 subgridMixingFactor = 100.
-# curvature acceleration coefficient
-# take curvature term into account in the gravity acceleration term
+# take curvature term into account in the gravity acceleration term for computing the friction force
+# (and gradient if curvAccInGradient=1)
 # 0 if deactivated, 1 if activated
-curvAcceleration = 0
+curvAccInFriction = 0
+# take curvature term into account in the tangential momentum equation
+# 0 if deactivated, 1 if activated
+curvAccInTangent = 0
+# if curvAccInFriction=1, take curvature term into account in the pressure gradient (if curvAccInGradient=1)
+curvAccInGradient = 0
 
 #++++++++++++Friction model
 # add the friction using an explicit formulation (1)
@@ -90,7 +96,10 @@ explicitFriction = 1
 frictModel = Coulomb
 #+++++++++++++General Friction parameters
 # tan of bed friction angle used for: samosAT, Coulomb, Voellmy
-mu = 0.45
+mu = 0.4
+#+++++++++++++Voellmy friction model
+xsi = 10000.
+
 
 [DFAPathGeneration_override]
 # use default DFAPathGeneration config as base configuration (True) and override following parameters
diff --git a/avaframe/ana1Tests/simiSolTest.py b/avaframe/ana1Tests/simiSolTest.py
index d85a8d085..e2a4f7edf 100644
--- a/avaframe/ana1Tests/simiSolTest.py
+++ b/avaframe/ana1Tests/simiSolTest.py
@@ -16,7 +16,6 @@
 import logging
 
 # local imports
-from avaframe.in3Utils import cfgUtils
 import avaframe.in3Utils.geoTrans as geoTrans
 import avaframe.com1DFA.com1DFA as com1DFA
 import avaframe.com1DFA.DFAtools as DFAtls
@@ -708,13 +707,20 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solSimi, fieldHeader, cfg
     vhErrorL2Array = np.zeros((len(fieldsList)))
     hErrorLMaxArray = np.zeros((len(fieldsList)))
     vhErrorLMaxArray = np.zeros((len(fieldsList)))
-    # get plots limits
-    limits = outAna1Plots.getPlotLimits(cfgSimi['SIMISOL'], fieldsList, fieldHeader)
     count = 0
+    tSave = cfgSimi['SIMISOL'].getfloat('tSave')
+    indT = min(np.searchsorted(timeList, tSave), min(len(timeList)-1, len(fieldsList)-1))
+    tSave = timeList[indT]
+    # get plots limits
+    if cfgSimi['SIMISOL'].getboolean('plotSequence'):
+        limits = outAna1Plots.getPlotLimits(cfgSimi['SIMISOL'], fieldsList, fieldHeader)
+    else:
+        limits = outAna1Plots.getPlotLimits(cfgSimi['SIMISOL'], [fieldsList[indT]], fieldHeader)
     # run the comparison routine for each saved time step
     for t, field in zip(timeList, fieldsList):
         # get similartiy solution h, u at required time step
         indTime = np.searchsorted(solSimi['time'], t)
+        log.debug('time simulation = %.2f, time simiSol = %.2f' % (t, solSimi['time'][indTime]))
         simiDict = getSimiSolParameters(solSimi, fieldHeader, indTime, cfgSimi['SIMISOL'], relTh, gravAcc)
         cellSize = fieldHeader['cellsize']
         cosAngle = simiDict['cos']
@@ -739,7 +745,7 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solSimi, fieldHeader, cfg
         log.debug("L2 %s error on the Flow Thickness at t=%.2f s is : %.4f" % (title, t, hErrorL2))
         log.debug("L2 %s error on the momentum at t=%.2f s is : %.4f" % (title, t, vhErrorL2))
         # Make all individual time step comparison plot
-        if cfgSimi['SIMISOL'].getboolean('plotSequence'):
+        if cfgSimi['SIMISOL'].getboolean('plotSequence') or t == tSave:
             outAna1Plots.saveSimiSolProfile(cfgMain, cfgSimi['SIMISOL'], field, limits, simiDict, t,
                                             fieldHeader, outDirTest, simHash)
             outAna1Plots.makeContourSimiPlot(avalancheDir, simHash, hNumerical, limits, simiDict, fieldHeader, t,
@@ -747,15 +753,13 @@ def analyzeResults(avalancheDir, fieldsList, timeList, solSimi, fieldHeader, cfg
         count = count + 1
 
     # Create result plots
-    tSave = cfgSimi['SIMISOL'].getfloat('tSave')
-    indT = min(np.searchsorted(timeList, tSave), min(len(timeList)-1, len(fieldsList)-1))
-    tSave = timeList[indT]
     indTime = np.searchsorted(solSimi['time'], tSave)
+    log.debug('time simulation = %.2f, time simiSol = %.2f' % (tSave, solSimi['time'][indTime]))
     simiDict = getSimiSolParameters(solSimi, fieldHeader, indTime, cfgSimi['SIMISOL'], relTh, gravAcc)
     outAna1Plots.plotSimiSolSummary(avalancheDir, timeList, fieldsList, fieldHeader, simiDict, hErrorL2Array,
                                     hErrorLMaxArray, vhErrorL2Array, vhErrorLMaxArray, outDirTest, simDFrow, simHash,
                                     cfgSimi['SIMISOL'])
-    if cfgSimi['SIMISOL'].getboolean('plotErrorTime') and len(timeList)>1:
+    if cfgSimi['SIMISOL'].getboolean('plotErrorTime') and len(timeList) > 1:
         outAna1Plots.plotErrorTime(timeList, hErrorL2Array, hErrorLMaxArray, vhErrorL2Array, vhErrorLMaxArray,
                                    cfgSimi['SIMISOL'].getboolean('relativError'), tSave, simHash, outDirTest)
 
@@ -773,7 +777,7 @@ def getReleaseThickness(avaDir, cfgSimi, demFile, csz):
     -----------
     avaDir: str
         path to avalanche directory
-    cfg: dict
+    cfgSimi: dict
         similarity solution confguration settings
     demFile: str
         path to DEM file
@@ -789,8 +793,7 @@ def getReleaseThickness(avaDir, cfgSimi, demFile, csz):
 
     # Read dem
     demOri = IOf.readRaster(demFile)
-    csz = cfg.getfloat('GENERAL', 'meshCellSize')
-    _, _, ncols, nrows = geoTrans.makeCoordGridFromHeader(demOri['header'], cellSizeNew=csz, larger=True)
+    _, _, ncols, nrows = geoTrans.makeCoordGridFromHeader(demOri['header'], cellSizeNew=csz, larger=False)
 
     xllc = demOri['header']['xllcenter']
     yllc = demOri['header']['yllcenter']
diff --git a/avaframe/ana1Tests/simiSolTestCfg.ini b/avaframe/ana1Tests/simiSolTestCfg.ini
new file mode 100644
index 000000000..5fa58b4f4
--- /dev/null
+++ b/avaframe/ana1Tests/simiSolTestCfg.ini
@@ -0,0 +1,132 @@
+### Config File - This file contains the main settings for the similarity solution test module
+## Set your parameters
+# This file will be overridden by local_simiSolTestCfg.ini if it exists
+# So copy this file to local_simiSolTestCfg.ini, adjust your variables there
+
+[SIMISOL]
+#sphKernel radius for the com1DFA simulations
+sphKernelRadius = 10|8|6|5|4|3
+
+# dimensioning parameters
+L_x = 80.
+L_y = 80.
+
+# release thickness
+relTh = 4.
+
+# bed friction angle
+bedFrictionAngle = 25.
+# internal friction angle
+internalFrictionAngle = 25.
+# plane inclination angle
+planeinclinationAngle = 35.
+
+# Flag earthpressure coefficients
+# if false takes 1
+flagEarth = False
+
+# save analysis plots at time step dtSol
+dtSol = 50.
+
+# which time step should be saved
+tSave = 20
+
+#++++Plotting parameters++++++
+# list of parameters to display in the summary plot (list of parameters separated by |)
+paramInfo = sphKernelRadius|aPPK|nPPK0|nPart
+# list of parameters to display in the convergence plot (list of parameters separated by |)
+paramConvInfo = sphKR0|nPPK0|cMax
+# plotting flags
+# only analyze and plot the tSave time step
+onlyLast = False
+# plot error function of time for each simulation
+plotErrorTime = False
+# plot individual figure for the h, hu and u error for each saved time step
+plotSequence = False
+# use relative difference
+relativError = True
+# when plotting, the domain extent is scaleCoef times bigger than the data extent
+scaleCoef = 1.02
+
+
+[com1DFA_override]
+# use default com1DFA config as base configuration (True) and override following parameters
+# if False and local is available use local
+defaultConfig = True
+#++++++++++++++++ Simulation type
+# list of simulations that shall be performed (null, ent, res, entres, available (use all available input data))
+simTypeList = null
+#+++++++++++++ Output++++++++++++
+# desired result Parameters (ppr, pfd, pfv, FD, FV, P, particles) - separated by |
+resType = FT|FV|Vx|Vy|Vz
+# saving time step, i.e.  time in seconds (first and last time step are always saved)
+# 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)
+# option 2: explicitly list all desired time steps (closest to actual computational time step) separated by | (example tSteps = 1|50.2|100)
+# NOTE: initial and last time step are always saved!
+tSteps = 0:1
+
+#+++++++++SNOW properties
+# True if release thickness should be read from shapefile file; if False - relTh read from ini file
+relThFromShp = False
+relThFromFile = True
+
+#++++++++++++Time stepping parameters
+# End time [s]
+tEnd = 30
+# to use a variable time step (time step depends on kernel radius)
+sphKernelRadiusTimeStepping = True
+# courant number
+cMax = 0.01
+
+#+++++++++++++SPH parameters
+# SPH gradient option
+# 0) No pressure gradients computed
+# 1) SamosAT style (no reprojecion on the surface, dz = 0 and gz is used)
+# 2) SamostAt but done in the cartesian coord system (will hopefully allow us to add the earth pressure coef)
+# 3) SamostAt but done in the local coord system (will hopefully allow us to add the earth pressure coef)
+# and this time reprojecion on the surface, dz not 0 and g3 is used
+sphOption = 2
+# sph kernel smoothing length [m]
+sphKernelRadius = 3
+# Choice of artificial viscosity
+# 0) No artificial viscosity
+# 1) SAMOS artificial viscosity
+# 2) Ata artificial viscosity
+viscOption = 1
+
+#++++++++++++++++ Particles
+# number of particles defined by: MPPDIR= mass per particle direct, MPPDH= mass per particles through release thickness,
+# MPPKR= mass per particles through number of particles per kernel radius
+massPerParticleDeterminationMethod = MPPKR
+# number of particles per kernel radius (if MPPKR is used)
+# is computed with: nPPK = nPPK0 * (sphKR/sphKR0)^aPPK
+# where sphKR is the sphKernelRadius specified further up
+# reference kernel radius
+sphKR0 = 5
+# reference number of particles per kernel radius
+nPPK0 = 15
+# variation of nppK exponent
+aPPK = 0|-0.5|-1|-2|-3
+
+# if splitOption=0
+# threshold for splitting particles, split if: mPart > massPerPart x thresholdMassSplit
+thresholdMassSplit = 5
+
+
+#+++++++++++++Mesh and interpolation
+# remesh the input DEM
+# expected mesh size [m]
+meshCellSize = 3
+
+#+++++++++++++Flow model parameters+++++++++
+# subgridMixingFactor
+subgridMixingFactor = 10
+
+#++++++++++++Friction model
+# friction type (samosAT, Coulomb)
+frictModel = Coulomb
+# add the friction using an explicit formulation (1)
+# 0 use an implicit
+explicitFriction = 1
+#+++++++++++++SamosAt friction model
+mu = 0.466307658154998593
diff --git a/avaframe/out3Plot/outAna1Plots.py b/avaframe/out3Plot/outAna1Plots.py
index e29c4daef..5f2f03463 100644
--- a/avaframe/out3Plot/outAna1Plots.py
+++ b/avaframe/out3Plot/outAna1Plots.py
@@ -57,7 +57,7 @@ def saveSimiSolProfile(cfgMain, cfgSimi, fields, limits, simiDict, tSave, header
     comSol['tSave'] = tSave
 
     # profile in flow direction
-    fig1 = plt.figure(figsize=(4*pU.figW, 2*pU.figH))
+    fig1 = plt.figure(figsize=(2.5*pU.figW, 1*pU.figH))
     fig1.suptitle('Similarty solution test, t = %.2f s (simulation %s)' % (tSave, simHash), fontsize=20)
     ax1 = plt.subplot2grid((1, 2), (0, 0))
     comSol['indFinal'] = int(nrows * 0.5) -1
@@ -76,7 +76,7 @@ def saveSimiSolProfile(cfgMain, cfgSimi, fields, limits, simiDict, tSave, header
 def makeContourSimiPlot(avalancheDir, simHash, fieldFT, limits, simiDict, fieldHeader, tSave, outDirTest):
     """
     """
-    fig, ax1 = plt.subplots(nrows=1, ncols=3, sharex=True, figsize=(pU.figW*4, pU.figH*2))
+    fig, ax1 = plt.subplots(nrows=1, ncols=3, sharex=True, figsize=(pU.figW*2, pU.figH))
     # make flow momentum comparison plot
     ax1 = plt.subplot2grid((1, 1), (0, 0))
     ax1 = addContour2Plot(ax1, fieldFT, simiDict, fieldHeader, limits, nLevels=16)
@@ -198,7 +198,7 @@ def _plotVariable(ax1, cfg, simiDict, comSol, limits, axis, particles=False):
 
 
 def plotSimiSolSummary(avalancheDir, timeList, fieldsList, fieldHeader, simiDict, hErrorL2Array, hErrorLMaxArray,
-                        vhErrorL2Array, vhErrorLMaxArray, outDirTest, simDFrow, simHash, cfgSimi):
+                       vhErrorL2Array, vhErrorLMaxArray, outDirTest, simDFrow, simHash, cfgSimi):
     """ Plot sumary figure of the similarity solution test
 
     Parameters
@@ -260,37 +260,36 @@ def plotSimiSolSummary(avalancheDir, timeList, fieldsList, fieldHeader, simiDict
     comSol['outDirTest'] = outDirTest
     comSol['tSave'] = tSave
 
-
     # create figures and plots
-    fig = plt.figure(figsize=(pU.figW*4, pU.figH*2))
+    fig = plt.figure(figsize=(pU.figW*2.5, pU.figH*2))
     fig.suptitle('Similarty solution test, t = %.2f s (simulation %s)' % (tSave, simHash), fontsize=30)
     # make comparison profile plot in flow direction
-    ax1 = plt.subplot2grid((2, 6), (0, 0), colspan=3)
-    comSol['indFinal'] = int(nrows * 0.5) -1
+    ax1 = plt.subplot2grid((2, 8), (0, 0), colspan=4)
+    comSol['indFinal'] = int(nrows * 0.5) - 1
     ax1, ax11 = _plotVariable(ax1, cfgSimi, simiDict, comSol, limits, 'xaxis')
     ax1.set_title('Profile in flow direction (y = 0 m)')
 
     # make flow momentum comparison plot
-    ax2 = plt.subplot2grid((2, 6), (0, 3), colspan=3)
+    ax2 = plt.subplot2grid((2, 8), (0, 4), colspan=4)
     comSol['indFinal'] = int(np.round((xCenter - xllc)/csz) + 1)
     ax2, ax22 = _plotVariable(ax2, cfgSimi, simiDict, comSol, limits, 'yaxis')
     ax2.set_title('Profile across flow direction (x = %.2f m)' % xCenter)
 
     # make bird view plot
-    ax3 = plt.subplot2grid((2, 6), (1, 0), colspan=2)
+    ax3 = plt.subplot2grid((2, 8), (1, 0), colspan=3)
     ax3 = addContour2Plot(ax3, fieldsList[indT]['FT'], simiDict, fieldHeader, limits)
     ax3.set_title(pU.cfgPlotUtils['nameFT'] + ' contours')
     pU.putAvaNameOnPlot(ax3, avalancheDir)
 
     # make error plot
-    ax4 = plt.subplot2grid((2, 6), (1, 2), colspan=2)
+    ax4 = plt.subplot2grid((2, 8), (1, 3), colspan=3)
     title = '\n between analytical solution and com1DFA'
     title = getTitleError(relativ, title)
     ax4.set_title(title)
     ax4, ax5 = addErrorTime(ax4, timeList, hErrorL2Array, hErrorLMaxArray, vhErrorL2Array, vhErrorLMaxArray,
                             relativ, tSave)
 
-    ax7 = plt.subplot2grid((2, 6), (1, 4), colspan=2)
+    ax7 = plt.subplot2grid((2, 8), (1, 6), colspan=2)
     ax7.axis("off")
     ax7.invert_yaxis()
     text = ''
@@ -408,7 +407,7 @@ def plotComparisonDam(cfg, simHash, fields0, fieldsT, header, solDam, tSave, lim
     # setup index for time of analyitcal solution
     indTime = np.searchsorted(solDam['tAna'], tSave)
 
-    fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, sharex=True, figsize=(pU.figW*4, pU.figH*2))
+    fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, sharex=True, figsize=(pU.figW*2, pU.figH*1))
     ax1 = _plotDamProfile(ax1, x, y, nx_loc, cfg, dataIniFT, dataAnaFT, solDam['xAna'], solDam['xMidAna'],
                                        solDam['hAna'], indTime, limits['maxFT'],  pU.cfgPlotUtils['nameFT'],
                                        pU.cfgPlotUtils['unitFT'])
@@ -490,11 +489,11 @@ def plotDamBreakSummary(avalancheDir, timeList, fieldsList, fieldHeader, solDam,
 
     """
     # Initialise DEM
-    demFile = gI.getDEMPath(avalancheDir)
-    demOri = IOf.readRaster(demFile, noDataToNan=True)
-    dem = com1DFA.initializeMesh(cfg['GENERAL'], demOri, cfg['GENERAL'].getint('methodMeshNormal'))
-    dem['header']['xllcenter'] = dem['originalHeader']['xllcenter']
-    dem['header']['yllcenter'] = dem['originalHeader']['yllcenter']
+    # demFile = gI.getDEMPath(avalancheDir)
+    # demOri = IOf.readRaster(demFile, noDataToNan=True)
+    # dem = com1DFA.initializeMesh(cfg['GENERAL'], demOri, cfg['GENERAL'].getint('methodMeshNormal'))
+    # dem['header']['xllcenter'] = dem['originalHeader']['xllcenter']
+    # dem['header']['yllcenter'] = dem['originalHeader']['yllcenter']
 
     cfgDam = cfg['DAMBREAK']
     phi = cfgDam.getfloat('phi')
@@ -548,7 +547,7 @@ def plotDamBreakSummary(avalancheDir, timeList, fieldsList, fieldHeader, solDam,
     limits = getPlotLimits(cfgDam, fieldsList[:indT], fieldHeader)
 
     # create figures and plots
-    fig = plt.figure(figsize=(pU.figW*4, pU.figH*2))
+    fig = plt.figure(figsize=(pU.figW*2.5, pU.figH*1.5))
     fig.suptitle('DamBreak test, t = %.2f s (simulation %s)' % (tSave, simHash), fontsize=30)
     # make flow thickness comparison plot
     ax1 = plt.subplot2grid((2, 6), (0, 0), colspan=2)
@@ -577,7 +576,7 @@ def plotDamBreakSummary(avalancheDir, timeList, fieldsList, fieldHeader, solDam,
     ax6 = plt.subplot2grid((2, 6), (1, 0), colspan=2)
     ax6, extent, cbar0, cs1 = outCom1DFA.addResult2Plot(ax6, fieldHeader, dataFT, 'FT')
     cbar0.ax.set_ylabel(pU.cfgPlotUtils['nameFT'])
-    ax6 = outCom1DFA.addDem2Plot(ax6, dem, what='slope', extent=extent)
+    # ax6 = outCom1DFA.addDem2Plot(ax6, dem, what='slope', extent=extent)
     rowsMin, rowsMax, colsMin, colsMax = pU.constrainPlotsToData(dataFT, fieldHeader['cellsize'],
                                                                  extentOption=True, constrainedData=False, buffer='')
     # draw rectangle corresponding to the error measurement domain
@@ -733,7 +732,7 @@ def addContour2Plot(ax1, fieldFT, simiDict, fieldHeader, limits, nLevels=9):
     # ax3.clabel(CS, inline=1, fontsize=8)
     h1,_ = cs1.legend_elements()
     h2,_ = cs2.legend_elements()
-    ax1.legend([h1[-1], h2[-1]], ['simulation', 'analytical'])
+    ax1.legend([h1[-1], h2[-1]], ['simulation', 'analytical'], loc='upper left')
     ax1.set_ylim([-widthY, widthY])
     ax1.set_xlim([-widthX+xCenter, widthX+xCenter])
     ax1.set_xlabel('x [m]')
@@ -741,7 +740,7 @@ def addContour2Plot(ax1, fieldFT, simiDict, fieldHeader, limits, nLevels=9):
     # add vertical and horizontal line showing the location of the profile plots cuts
     ax1.axvline(xCenter, color='b', linestyle=':')
     ax1.axhline(0, color='r', linestyle=':')
-    ax1.text(xCenter+5, -widthY, "x = %.2f m" % (xCenter), color='b')
+    ax1.text(xCenter+5, widthY, "x = %.2f m" % (xCenter), color='b', verticalalignment='top')
     ax1.text(-widthX+xCenter, 5, "y = 0 m", color='r')
 
     return ax1
@@ -774,11 +773,12 @@ def plotErrorConvergence(simDF, outDirTest, cfgSimi, xField, yField, coloredBy,
     fit: boolean
         if True add power law regression
     """
+    paramConvInfo = cfgSimi['paramConvInfo'].split('|')
     tSave = simDF['timeError'][0]
     relativ = cfgSimi.getboolean('relativError')
     cmap, _, ticks, norm = pU.makeColorMap(pU.cmapAvaframeCont, min(simDF[coloredBy]), max(simDF[coloredBy]),
                                            continuous=pU.contCmap)
-    fig1, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(3*pU.figW, 2*pU.figH))
+    fig1, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(2*pU.figW, pU.figH))
     # get the sizing function
     sizeList = simDF[sizedBy].unique()
     lenSize = len(sizeList)
@@ -808,20 +808,25 @@ def plotErrorConvergence(simDF, outDirTest, cfgSimi, xField, yField, coloredBy,
             ax1.plot(xArray, p0H*xArray**p1H, color=color)
             if np.size(rSquaredH) == 0:
                 rSquaredH = np.nan
-            log.debug('power law fit sphKernelRadius = %.2f m: hErrorL2 = %.1f * Npart^{%.2f}, r=%.2f' %
-                     (colorValue, p0H, p1H, rSquaredH))
+            log.debug('power law fit sphKernelRadius = %.2f m: %s = %.1f * Npart^{%.2f}, r=%.2f' %
+                     (colorValue, yField, p0H, p1H, rSquaredH))
 
     if logScale:
         ax1.set_yscale('log')
         ax1.set_xscale('log')
 
     fig1.suptitle('Convergence of DFA simulation for the similarity solution test at t = %.2fs' % tSave)
-    ax1.set_title(getTitleError(relativ, r' L2 on flow thickness'))
+    if yField == 'vhErrorL2':
+        what = getLabel(' on flow momentum ', '', dir='', vert=True)
+    else:
+        what = r' on flow thickness $h$'
+    ax1.set_title(getTitleError(relativ, what))
     ax1.set_xlabel(xField)
-    ax1.set_ylabel(getTitleError(relativ, r' L2 on flow thickness'))
+    ax1.set_ylabel(getTitleError(relativ, what))
     if lenColor<=10:
         lenColor = None
-    legend1 = ax1.legend(*scatter1.legend_elements(num=lenColor), loc="upper center", title=coloredBy)
+    legend1 = ax1.legend(*scatter1.legend_elements(num=lenColor), bbox_to_anchor=(1, 0.95), loc="upper right",
+                         title=coloredBy)
     ax1.add_artist(legend1)
 
     # produce a legend with a cross section of sizes from the scatter
@@ -829,9 +834,16 @@ def plotErrorConvergence(simDF, outDirTest, cfgSimi, xField, yField, coloredBy,
         lenSize = None
     kw = dict(prop="sizes", color=scatter1.cmap(0.7),
           func=lambda s: (s-10)*(maxSize - minSize)/70 + minSize)
-    legend3 = ax1.legend(*scatter1.legend_elements(num=lenSize, **kw), loc="upper right", title=sizedBy)
-    ax1.grid(color='grey', linestyle='-', linewidth=0.25, alpha=0.5)
-    ax1.grid(color='grey', which='minor', linestyle='--', linewidth=0.25, alpha=0.5)
+    legend3 = ax1.legend(*scatter1.legend_elements(num=lenSize, **kw), bbox_to_anchor=(0.8, 0.95), loc="upper right",
+                         title=sizedBy)
+    ax1.grid(color='grey', linestyle='-', linewidth=0.5, alpha=0.75)
+    ax1.grid(color='grey', which='minor', linestyle='--', linewidth=0.5, alpha=0.5)
+
+    text = ''
+    for param in paramConvInfo:
+        text = text + (param + ' = %.2f' % simDF[param][0]) + '\n'
+    ax1.text(0.1, 0.1, text, transform=ax1.transAxes, ha='left', va='bottom', fontsize=pU.fs)
+
     pU.saveAndOrPlot({'pathResult': outDirTest}, 'ErrorLog%ds' % int(tSave), fig1)
     return fig1, ax1
 
@@ -909,9 +921,9 @@ def plotPresentation(simDF, outDirTest, cfgSimi, xField, yField, coloredBy, size
         scatter1 = ax1.scatter(simDFNew[xField], simDFNew[yField], c=simDFNew[coloredBy], s=sizeList,
                                cmap=cmap, norm=norm, marker=pU.markers[0], alpha=1)
         fig1.suptitle('Convergence of DFA simulation for the similarity solution test at t = %.2fs' % tSave)
-        ax1.set_title(getTitleError(relativ, r' L2 on flow thickness'))
+        ax1.set_title(getTitleError(relativ, r' on flow thickness'))
         ax1.set_xlabel(xField)
-        ax1.set_ylabel(getTitleError(relativ, r' L2 on flow thickness'))
+        ax1.set_ylabel(getTitleError(relativ, r' on flow thickness'))
         legend1 = ax1.legend(*scatter1.legend_elements(num=lenColor), loc="upper center", title=coloredBy)
         ax1.add_artist(legend1)
 
@@ -1088,9 +1100,9 @@ def getPlotLimits(cfgSimi, fieldsList, fieldHeader):
 def getTitleError(relativ, ending=''):
     """Get error plot title (relativ error or not?)"""
     if relativ:
-        return 'Relative error difference' + ending
+        return 'Relative L2 deviation' + ending
     else:
-        return 'Error difference' + ending
+        return 'L2 deviation' + ending
 
 
 def getLabel(start, end, dir='', vert=True):
diff --git a/avaframe/out3Plot/plotUtils.py b/avaframe/out3Plot/plotUtils.py
index cebcce02a..37a226938 100644
--- a/avaframe/out3Plot/plotUtils.py
+++ b/avaframe/out3Plot/plotUtils.py
@@ -68,6 +68,12 @@
 matplotlib.rcParams['xtick.major.width'] = 1.0
 matplotlib.rcParams['xtick.major.size'] = 3
 matplotlib.rcParams['xtick.labelsize'] = cfg['tickLabelSize']
+matplotlib.rcParams['xtick.direction'] = 'in'
+matplotlib.rcParams['xtick.top'] = True
+matplotlib.rcParams['xtick.bottom'] = True
+matplotlib.rcParams['ytick.direction'] = 'in'
+matplotlib.rcParams['ytick.left'] = True
+matplotlib.rcParams['ytick.right'] = True
 matplotlib.rcParams['ytick.color'] = 'grey'
 matplotlib.rcParams['ytick.major.width'] = 1.0
 matplotlib.rcParams['ytick.major.size'] = 3
diff --git a/avaframe/runScripts/runAnalyzeDamBreak.py b/avaframe/runScripts/runAnalyzeDamBreak.py
deleted file mode 100644
index 0d630a550..000000000
--- a/avaframe/runScripts/runAnalyzeDamBreak.py
+++ /dev/null
@@ -1,69 +0,0 @@
-"""
-    Run script for running the dam break problem on a sloping bed and compare to simulation results
-"""
-
-# Load modules
-import numpy as np
-import pathlib
-import pandas as pd
-
-# Local imports
-from avaframe.com1DFA import com1DFA
-from avaframe.in3Utils import fileHandlerUtils as fU
-from avaframe.ana1Tests import damBreak
-from avaframe.in3Utils import cfgUtils
-from avaframe.in3Utils import logUtils
-import avaframe.out3Plot.outAna1Plots as outAna1Plots
-
-
-
-# log file name; leave empty to use default runLog.log
-logName = 'runAnalyzeDamBreakProblem'
-
-# Load settings from general configuration file
-cfgMain = cfgUtils.getGeneralConfig()
-avalancheDir = 'data/avaDamBreak'
-cfgMain['MAIN']['avalancheDir'] = avalancheDir
-
-# Start logging
-log = logUtils.initiateLogger(avalancheDir, logName)
-log.info('MAIN SCRIPT')
-log.info('Current avalanche: %s', avalancheDir)
-
-# create output directory for test result plots
-outDirTest = pathlib.Path(avalancheDir, 'Outputs', 'ana1Tests')
-fU.makeADir(outDirTest)
-
-# Load configuration
-damBreakCfg = pathlib.Path(avalancheDir, 'Inputs', 'damBreak_com1DFACfg.ini')
-cfg = cfgUtils.getModuleConfig(com1DFA, damBreakCfg)
-cfgGen = cfg['GENERAL']
-
-# Load flow thickness from analytical solution
-solDam = damBreak.damBreakSol(avalancheDir, cfgMain, cfg, outDirTest)
-
-# create dataFrame of results
-simDF, _ = cfgUtils.readAllConfigurationInfo(avalancheDir)
-
-
-# if the analysis already exists and you only want to replot uncomment this (and put the good result name)
-# pathToResults = pathlib.Path(avalancheDir, 'Outputs', 'ana1Tests', 'results5.p')
-# if pathToResults.is_file():
-#     simDF = pd.read_pickle(pathToResults)
-simDF = damBreak.postProcessDamBreak(avalancheDir, cfgMain, cfg, simDF, solDam, outDirTest)
-
-# make convergence plot (if you add the fiting lines, make sure only the coloredBy and sizedBy parameters are varied)
-fig1, ax1 = outAna1Plots.plotErrorConvergence(simDF, outDirTest, cfg['DAMBREAK'], 'nPart', 'vhErrorL2',
-                          'aPPK', 'nPPK0', logScale=True, fit=False)
-
-# make convergence plot
-outAna1Plots.plotTimeCPULog(simDF, outDirTest, cfg['DAMBREAK'], 'nPart', 'aPPK', 'nPPK0')
-
-# do some filtering for the presentation plot
-simDF = simDF[simDF['sphKernelRadius']==3]
-
-# make convergence plot (if you add the fiting lines, make sure only the coloredBy and sizedBy parameters are varied)
-# same as plotErrorConvergence but adds the points corresponding to different coloredBy values one after the others
-# and saves itermediate plots
-fig1, ax1 = outAna1Plots.plotPresentation(simDF, outDirTest, cfg['DAMBREAK'], 'nPart', 'hErrorL2',
-                          'aPPK', 'nPPK0', logScale=True, fit=True)
diff --git a/avaframe/runScripts/runAnalyzeSimilaritySol.py b/avaframe/runScripts/runAnalyzeSimilaritySol.py
index 63c4aee4d..a2e7204ee 100644
--- a/avaframe/runScripts/runAnalyzeSimilaritySol.py
+++ b/avaframe/runScripts/runAnalyzeSimilaritySol.py
@@ -11,7 +11,6 @@
 """
 
 import pathlib
-import pandas as pd
 
 # Local imports
 import avaframe.in3Utils.initializeProject as initProj
@@ -25,7 +24,7 @@
 
 # +++++++++SETUP CONFIGURATION++++++++++++++++++++++++
 # log file name; leave empty to use default runLog.log
-logName = 'runSimilarityTest'
+logName = 'runAnalyzeSimilarityTest'
 
 # Load general configuration
 cfgMain = cfgUtils.getGeneralConfig()
@@ -72,7 +71,7 @@
 
 # make convergence plot (if you add the fiting lines, make sure only the coloredBy and sizedBy parameters are varied)
 fig1, ax1 = outAna1Plots.plotErrorConvergence(simDF, outDirTest, simiSolCfg['SIMISOL'], 'nPart', 'vhErrorL2',
-                                              'aPPK', 'nPPK0', logScale=True, fit=False)
+                                              'aPPK', 'sphKernelRadius', logScale=True, fit=True)
 
 # make convergence plot
 outAna1Plots.plotTimeCPULog(simDF, outDirTest, simiSolCfg['SIMISOL'], 'nPart', 'sphKernelRadius', 'nPPK0')
@@ -84,6 +83,6 @@
 # same as plotErrorConvergence but adds the points corresponding to different coloredBy values one after the others
 # and saves itermediate plots
 fig1, ax1 = outAna1Plots.plotPresentation(simDF, outDirTest, simiSolCfg['SIMISOL'], 'nPart', 'hErrorL2',
-                                          'aPPK', 'nPPK0', logScale=True, fit=True)
+                                          'aPPK', 'sphKernelRadius', logScale=True, fit=True)
 
 outAna1Plots.plotTimeCPULog(simDF, outDirTest, simiSolCfg['SIMISOL'], 'nPart', 'aPPK', 'nPPK0')
diff --git a/avaframe/runScripts/runDamBreak.py b/avaframe/runScripts/runDamBreak.py
index a6672208a..568abeddc 100644
--- a/avaframe/runScripts/runDamBreak.py
+++ b/avaframe/runScripts/runDamBreak.py
@@ -1,5 +1,5 @@
 """
-    Run script for running the dam break problem on a sloping bed and compare to simulation results
+    Run script for running the dam break test DFA simulations on an inclined plane and compare to simulation results
 """
 
 # Load modules
@@ -7,14 +7,23 @@
 from configupdater import ConfigUpdater
 
 # Local imports
+import avaframe.in3Utils.fileHandlerUtils as fU
 from avaframe.com1DFA import com1DFA
 import avaframe.in3Utils.initializeProject as initProj
 from avaframe.in3Utils import cfgUtils
+from avaframe.in3Utils import cfgHandling
 from avaframe.in3Utils import logUtils
+import avaframe.ana1Tests.damBreakTest as damBreakTest
+import avaframe.out3Plot.outAna1Plots as outAna1Plots
 
 
+# +++++++++ REQUIRED ++++++++++++++++++++++++
 # log file name; leave empty to use default runLog.log
-logName = 'runDamBreakProblem'
+logName = 'runDamBreakTest'
+# if left empty, use the damBreakTestCfg.ini and local_damBreakTestCfg.ini configuration files
+# use 'ana1Tests/figConvergence_damBreakTestCfg.ini' to produce the similarity solution plots from the Theory Paper
+fileOverride = ''
+# ++++++++++++++++++++++++++++++
 
 # Load settings from general configuration file
 cfgMain = cfgUtils.getGeneralConfig()
@@ -22,23 +31,69 @@
 cfgMain['MAIN']['avalancheDir'] = avalancheDir
 
 # Clean input directory(ies) of old work and output files
-# initProj.cleanModuleFiles(avalancheDir, com1DFA)
+initProj.cleanSingleAvaDir(avalancheDir, keep=logName, deleteOutput=False)
+# setup work folder
+workPath = pathlib.Path(avalancheDir, 'Work', 'ana1Tests', 'damBreakTest')
+fU.makeADir(workPath)
+# create output directory for test result plots
+outDirTest = pathlib.Path(avalancheDir, 'Outputs', 'ana1Tests', 'damBreakTest')
+fU.makeADir(outDirTest)
 
 # Start logging
 log = logUtils.initiateLogger(avalancheDir, logName)
 log.info('MAIN SCRIPT')
 log.info('Current avalanche: %s', avalancheDir)
 
-# Load configuration
-damBreakCfg = pathlib.Path(avalancheDir, 'Inputs', 'damBreak_com1DFACfg.ini')
+# Load configuration for dam break test
+damBreakCfg = cfgUtils.getModuleConfig(damBreakTest, fileOverride=fileOverride)
+# ++++++++++ set configurations for all the used modules and override ++++++++++++
+# get comDFA configuration and save to file
+com1DFACfg = cfgUtils.getModuleConfig(com1DFA, fileOverride='', modInfo=False, toPrint=False,
+                                      onlyDefault=damBreakCfg['com1DFA_override']['defaultConfig'])
+com1DFACfg, damBreakCfg = cfgHandling.applyCfgOverride(com1DFACfg, damBreakCfg, com1DFA, addModValues=False)
+com1DFACfgFile = cfgUtils.writeCfgFile(avalancheDir, com1DFA, com1DFACfg, fileName='com1DFA_settings',
+                                       filePath=workPath)
 
-# Overwriting the sphKernelRadius and meshCellSize from the ini file with the values specified in []
-for sphKernelRadius in [5, 3]:
+# run DFA simulations
+sphKernelRadiusList = damBreakCfg['DAMBREAK']['sphKernelRadius'].split('|')
+sphKernelRadiusList = [float(i) for i in sphKernelRadiusList]
+for sphKernelRadius in sphKernelRadiusList:
     updater = ConfigUpdater()
-    updater.read(damBreakCfg)
+    updater.read(com1DFACfgFile)
     updater['GENERAL']['sphKernelRadius'].value = sphKernelRadius
     updater['GENERAL']['meshCellSize'].value = sphKernelRadius
     updater.update_file()
 
     # call com1DFA to perform simulation - provide configuration file and release thickness function
-    dem, plotDict, reportDictList, simDF = com1DFA.com1DFAMain(avalancheDir, cfgMain, cfgFile=damBreakCfg)
+    initProj.cleanModuleFiles(avalancheDir, com1DFA, deleteOutput=False)
+    dem, plotDict, reportDictList, simDF = com1DFA.com1DFAMain(avalancheDir, cfgMain, cfgFile=com1DFACfgFile)
+
+# analyze simulations
+# Load configuration info of all com1DFA simulations
+simDF, _ = cfgUtils.readAllConfigurationInfo(avalancheDir)
+
+# Load flow thickness from analytical solution
+solDam = damBreakTest.damBreakSol(avalancheDir, cfgMain, damBreakCfg['DAMBREAK'], com1DFACfg, outDirTest)
+
+
+# if the analysis already exists and you only want to replot uncomment this (and put the good result name)
+# pathToResults = pathlib.Path(avalancheDir, 'Outputs', 'ana1Tests', 'damBreakTest', 'results5.p')
+# if pathToResults.is_file():
+#     simDF = pd.read_pickle(pathToResults)
+simDF = damBreakTest.postProcessDamBreak(avalancheDir, cfgMain, damBreakCfg, simDF, solDam, outDirTest)
+
+# do some filtering for the presentation plot
+# simDF = simDF[simDF['sphKernelRadius']==3]
+
+# make convergence plot (if you add the fiting lines, make sure only the coloredBy and sizedBy parameters are varied)
+fig1, ax1 = outAna1Plots.plotErrorConvergence(simDF, outDirTest, damBreakCfg['DAMBREAK'], 'nPart', 'vhErrorL2',
+                                              'aPPK', 'sphKernelRadius', logScale=True, fit=True)
+
+# # make convergence plot
+# outAna1Plots.plotTimeCPULog(simDF, outDirTest, damBreakCfg['DAMBREAK'], 'nPart', 'aPPK', 'sphKernelRadius')
+#
+# # make convergence plot (if you add the fiting lines, make sure only the coloredBy and sizedBy parameters are varied)
+# # same as plotErrorConvergence but adds the points corresponding to different coloredBy values one after the others
+# # and saves itermediate plots
+# fig1, ax1 = outAna1Plots.plotPresentation(simDF, outDirTest, damBreakCfg['DAMBREAK'], 'nPart', 'vhErrorL2',
+#                                           'aPPK', 'sphKernelRadius', logScale=True, fit=True)
diff --git a/avaframe/runScripts/runEnergyLineTest.py b/avaframe/runScripts/runEnergyLineTest.py
index 83d7ada33..501a1b7d9 100644
--- a/avaframe/runScripts/runEnergyLineTest.py
+++ b/avaframe/runScripts/runEnergyLineTest.py
@@ -22,6 +22,10 @@
 # +++++++++REQUIRED+++++++++++++
 # log file name; leave empty to use default runLog.log
 logName = 'runEnergyLineTest'
+# if left empty, use the energyLineTestCfg.ini and local_energyLineTestCfg.ini configuration files
+# use 'ana1Tests/figAnalytic_energyLineTestCfg.ini' or 'ana1Tests/figCurvature_energyLineTestCfg.ini'
+# or 'ana1Tests/figPressure_energyLineTestCfg.ini' to produce the energy line plots from the Theory Paper
+fileOverride = ''
 # ++++++++++++++++++++++++++++++
 
 # Load avalanche directory from general configuration file
@@ -38,7 +42,7 @@
 iP.cleanSingleAvaDir(avalancheDir, keep=logName, deleteOutput=False)
 workPath = pathlib.Path(avalancheDir, 'Work', 'energyLineTest')
 fU.makeADir(workPath)
-energyLineTestCfg = cfgUtils.getModuleConfig(energyLineTest)
+energyLineTestCfg = cfgUtils.getModuleConfig(energyLineTest, fileOverride=fileOverride)
 
 # ++++++++++ set configurations for all the used modules and override ++++++++++++
 # get comDFA configuration and save to file
diff --git a/avaframe/runScripts/runRotationTest.py b/avaframe/runScripts/runRotationTest.py
index 0e2695439..6cc5347d9 100644
--- a/avaframe/runScripts/runRotationTest.py
+++ b/avaframe/runScripts/runRotationTest.py
@@ -23,6 +23,10 @@
 # +++++++++REQUIRED+++++++++++++
 # log file name; leave empty to use default runLog.log
 logName = 'runRotationTest'
+# if left empty, use the rotationTestCfg.ini and local_rotationTestCfg.ini configuration files
+# use 'ana1Tests/figPressure_rotationTestCfg.ini' or 'ana1Tests/figAnalytic_rotationTestCfg.ini' to produce
+# the rotation test plots from the Theory Paper
+fileOverride = ''
 # ++++++++++++++++++++++++++++++
 
 # Load avalanche directory from general configuration file
@@ -40,7 +44,7 @@
 iP.cleanSingleAvaDir(avalancheDir, keep=logName, deleteOutput=False)
 workPath = pathlib.Path(avalancheDir, 'Work', 'energyLineTest')
 fU.makeADir(workPath)
-rotationTestCfg = cfgUtils.getModuleConfig(rotationTest)
+rotationTestCfg = cfgUtils.getModuleConfig(rotationTest, fileOverride=fileOverride)
 
 # ++++++++++ set configurations for all the used modules and override ++++++++++++
 # get comDFA configuration and save to file
diff --git a/avaframe/runScripts/runSimilaritySol.py b/avaframe/runScripts/runSimilaritySol.py
index cf811ea69..27f143765 100644
--- a/avaframe/runScripts/runSimilaritySol.py
+++ b/avaframe/runScripts/runSimilaritySol.py
@@ -1,5 +1,5 @@
 """
-    Run com1DFA kernel and compare it tothe similarity solution
+    Run com1DFA kernel and compare it to the similarity solution
     This script computes the similarity solution for a gliding avalanche on
     a inclined plane according to similarity solution from :
     Hutter, K., Siegel, M., Savage, S.B. et al.
@@ -10,6 +10,7 @@
 """
 
 import pathlib
+import pandas as pd
 from configupdater import ConfigUpdater
 
 # Local imports
@@ -21,12 +22,18 @@
 import avaframe.com1DFA.com1DFA as com1DFA
 from avaframe.in1Data import getInput as gI
 import avaframe.ana1Tests.simiSolTest as simiSolTest
+import avaframe.out3Plot.outAna1Plots as outAna1Plots
 
 
-# +++++++++SETUP CONFIGURATION++++++++++++++++++++++++
+# +++++++++ REQUIRED ++++++++++++++++++++++++
 # log file name; leave empty to use default runLog.log
 logName = 'runSimilarityTest'
 
+# if left empty, use the simiSolTestCfg.ini and local_simiSolTestCfg.ini configuration files
+# use 'ana1Tests/figConvergence_simiSolTestCfg.ini' to produce the similarity solution plots from the Theory Paper
+fileOverride = ''
+# ++++++++++++++++++++++++++++++
+
 # Load general configuration
 cfgMain = cfgUtils.getGeneralConfig()
 avalancheDir = 'data/avaSimilaritySol'
@@ -46,7 +53,7 @@
 log.info('Current avalanche: %s', avalancheDir)
 
 # Load configuration for similarity solution test
-simiSolCfg = cfgUtils.getModuleConfig(simiSolTest)
+simiSolCfg = cfgUtils.getModuleConfig(simiSolTest, fileOverride=fileOverride)
 # ++++++++++ set configurations for all the used modules and override ++++++++++++
 # get comDFA configuration and save to file
 com1DFACfg = cfgUtils.getModuleConfig(com1DFA, fileOverride='', modInfo=False, toPrint=False,
@@ -72,3 +79,34 @@
     # (may be multiple sims)
     initProj.cleanModuleFiles(avalancheDir, com1DFA, deleteOutput=False)
     _, _, _, simDF = com1DFA.com1DFAMain(avalancheDir, cfgMain, cfgFile=com1DFACfgFile)
+
+
+# analyze simulations
+# Load configuration info of all com1DFA simulations
+simDF, _ = cfgUtils.readAllConfigurationInfo(avalancheDir)
+
+# compute the similartiy solution (this corresponds to our reference)
+log.info('Computing similarity solution')
+solSimi = simiSolTest.mainSimilaritySol(simiSolCfg['SIMISOL'], com1DFACfg)
+
+# if the analysis already exists and you only want to replot uncomment this (and put the good result name)
+# pathToResults = pathlib.Path(avalancheDir, 'Outputs', 'ana1Tests', 'simiSolTest', 'results20.p')
+# if pathToResults.is_file():
+#     simDF = pd.read_pickle(pathToResults)
+simDF = simiSolTest.postProcessSimiSol(avalancheDir, cfgMain, simiSolCfg, simDF, solSimi, outDirTest)
+
+# do some filtering before plotting
+# simDF = simDF[simDF['nPPK0'] == 30]
+
+# make convergence plot (if you add the fiting lines, make sure only the coloredBy and sizedBy parameters are varied)
+fig1, ax1 = outAna1Plots.plotErrorConvergence(simDF, outDirTest, simiSolCfg['SIMISOL'], 'nPart', 'vhErrorL2',
+                                              'aPPK', 'sphKernelRadius', logScale=True, fit=True)
+
+# # make cpu time plot
+# outAna1Plots.plotTimeCPULog(simDF, outDirTest, simiSolCfg['SIMISOL'], 'nPart', 'sphKernelRadius', 'sphKernelRadius')
+#
+# # make convergence plot (if you add the fiting lines, make sure only the coloredBy and sizedBy parameters are varied)
+# # same as plotErrorConvergence but adds the points corresponding to different coloredBy values one after the others
+# # and saves itermediate plots
+# fig1, ax1 = outAna1Plots.plotPresentation(simDF, outDirTest, simiSolCfg['SIMISOL'], 'nPart', 'hErrorL2',
+#                                           'aPPK', 'sphKernelRadius', logScale=True, fit=True)