diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx
index 7e8b8698d..f6fc61ff5 100644
--- a/avaframe/com1DFA/DFAfunctionsCython.pyx
+++ b/avaframe/com1DFA/DFAfunctionsCython.pyx
@@ -194,7 +194,6 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
       indCellY = indYDEM[k]
       # deduce area
       areaPart = m / (h * rho)
-
       # get cell and weights
       Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
 
@@ -2293,3 +2292,75 @@ def computeIniMovement(cfg, particles, dem, dT, fields):
   particles['m'] = np.asarray(mass)
 
   return particles, force
+
+
+def updateInitialVelocity(cfg, particles, dem, float velocityMag):
+  """
+  update particles' velocity in direction of the steepest descent
+
+  Parameters
+  ------------
+  cfg: configparser
+    configuration for DFA simulation
+  particles : dict
+    particles dictionary at t
+  dem : dict
+    dictionary with dem information
+  velocityMag: float
+    velocity of the particles
+
+  Returns
+  ------------
+  particles : dict
+    particles dictionary at t with updated velocity
+  """
+
+  cdef double[:] xArray = particles['x']
+  cdef double[:] yArray = particles['y']
+  cdef double[:] uxArray = particles['ux']
+  cdef double[:] uyArray = particles['uy']
+  cdef double[:] uzArray = particles['uz']
+  cdef double[:] velocityMagArray = particles['velocityMag']
+  cdef int nPart = particles['nPart']
+  cdef int nrows = dem['header']['nrows']
+  cdef int ncols = dem['header']['ncols']
+  cdef double[:, :] nxArray = dem['Nx']
+  cdef double[:, :] nyArray = dem['Ny']
+  cdef double[:, :] nzArray = dem['Nz']
+  cdef int interpOption = cfg.getint('interpOption')
+  cdef double csz = dem['header']['cellsize']
+
+  cdef double x, y, z
+  cdef double nx, ny, nz
+  cdef double tangentialTopoX, tangentialTopoY, tangentialTopoZ, tangentialTopoXNorm, tangentialTopoYNorm, tangentialTopoZNorm
+  cdef int Lx0, Ly0, iCell
+  cdef double w[4]
+  cdef int k
+
+  for k in range(nPart):
+    x = xArray[k]
+    y = yArray[k]
+
+    # get cell and weights
+    Lx0, Ly0, iCell, w[0], w[1], w[2], w[3] = DFAtlsC.getCellAndWeights(x, y, ncols, nrows, csz, interpOption)
+
+    # get normal at the particle location
+    nx, ny, nz = DFAtlsC.getVector(Lx0, Ly0, w[0], w[1], w[2], w[3], nxArray, nyArray, nzArray)
+    nx, ny, nz = DFAtlsC.normalize(nx, ny, nz)
+    # get vector in the tangent plane that points in the direction gravity would pull a particle along the surface
+    tangentialTopoX = nx * nz
+    tangentialTopoY = ny * nz
+    tangentialTopoZ = nz * nz -1
+    tangentialTopoXNorm, tangentialTopoYNorm, tangentialTopoZNorm = DFAtlsC.normalize(tangentialTopoX, tangentialTopoY, tangentialTopoZ)
+
+    uxArray[k] = tangentialTopoXNorm * velocityMag
+    uyArray[k] = tangentialTopoYNorm * velocityMag
+    uzArray[k] = tangentialTopoZNorm * velocityMag
+
+    velocityMagArray[k] = DFAtlsC.norm(uxArray[k], uyArray[k], uzArray[k])
+  particles['ux'] = np.asarray(uxArray)
+  particles['uy'] = np.asarray(uyArray)
+  particles['uz'] = np.asarray(uzArray)
+  particles['velocityMag'] = np.asarray(velocityMagArray)
+
+  return particles
\ No newline at end of file
diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py
index 55b9740a6..1d01e5406 100644
--- a/avaframe/com1DFA/com1DFA.py
+++ b/avaframe/com1DFA/com1DFA.py
@@ -53,6 +53,7 @@
 from avaframe.com1DFA import checkCfg
 from avaframe.ana5Utils import distanceTimeAnalysis as dtAna
 import avaframe.out3Plot.outDistanceTimeAnalysis as dtAnaPlots
+import avaframe.com1DFA.debrisFunctions as debF
 import threading
 
 #######################################
@@ -562,6 +563,8 @@ def prepareInputData(inputSimFiles, cfg):
         - secondaryRelFile : str, path to secondaryRelease file
         - entFiles : str, path to entrainment file
         - resFile : str, path to resistance file
+        - hydrographFile: str, path to hydrograph polygon file
+        - hydrographCsv: str, path to hydrograph values (csv-)file
         - entResInfo : flag dict
         flag if Yes entrainment and/or resistance areas found and used for simulation
         flag True if a Secondary Release file found and activated
@@ -585,6 +588,7 @@ def prepareInputData(inputSimFiles, cfg):
         - resLine : dict, resistance line dictionary
         - entrainmentArea : str, entrainment file name
         - resistanceArea : str, resistance file name
+        - hydrographAreaLine: dict, hydrograph line dictionary
         - entResInfo : flag dict
         flag if Yes entrainment and/or resistance areas found and used for simulation
         flag True if a Secondary Release file found and activated
@@ -741,6 +745,11 @@ def prepareInputData(inputSimFiles, cfg):
     else:
         damLine = None
 
+    if cfg["GENERAL"].getboolean("hydrograph"):
+        hydrLine = debF.preparehydrographAreaLine(inputSimFiles, demOri, cfg)
+    else:
+        hydrLine = None
+
     inputSimLines = {
         "releaseLine": releaseLine,
         "secondaryReleaseLine": secondaryReleaseLine,
@@ -756,6 +765,7 @@ def prepareInputData(inputSimFiles, cfg):
         "xiFile": inputSimFiles["xiFile"],
         "kFile": inputSimFiles["kFile"],
         "tauCFile": inputSimFiles["tauCFile"],
+        "hydrographAreaLine": hydrLine,
     }
 
     return demOri, inputSimLines
@@ -1613,6 +1623,8 @@ def getRelThFromPart(cfg, releaseLine, relThField, thName):
 
     if len(relThField) != 0:
         relThForPart = np.amax(relThField)
+    elif releaseLine["type"] == "Hydrograph":
+        relThForPart = releaseLine["thickness"]
     elif cfg.getboolean("%sThFromFile" % thName):
         relThForPart = np.amax(np.asarray(releaseLine["thickness"], dtype=float))
     else:
@@ -1654,26 +1666,26 @@ def initializeFields(cfg, dem, particles, releaseLine):
     nrows = header["nrows"]
     # initialize fields
     fields = {}
-    fields["pfv"] = np.zeros((nrows, ncols))
-    fields["pft"] = np.zeros((nrows, ncols))
-    fields["FV"] = np.zeros((nrows, ncols))
-    fields["FT"] = np.zeros((nrows, ncols))
-    fields["FM"] = np.zeros((nrows, ncols))
-    fields["Vx"] = np.zeros((nrows, ncols))
-    fields["Vy"] = np.zeros((nrows, ncols))
-    fields["Vz"] = np.zeros((nrows, ncols))
-    fields["dmDet"] = np.zeros((nrows, ncols))
-    fields["FTStop"] = np.zeros((nrows, ncols))
-    fields["FTDet"] = np.zeros((nrows, ncols))
-    fields["FTEnt"] = np.zeros((nrows, ncols))
-    fields["sfcChange"] = np.zeros((nrows, ncols))
-    fields["sfcChangeTotal"] = np.zeros((nrows, ncols))
-    fields["demAdapted"] = np.zeros((nrows, ncols))
+    fields["pfv"] = np.zeros((nrows, ncols))  # peak flow velocity [m/s]
+    fields["pft"] = np.zeros((nrows, ncols))  # peal flow thickness [m]
+    fields["FV"] = np.zeros((nrows, ncols))  # flow velocity [m/s]
+    fields["FT"] = np.zeros((nrows, ncols))  # flow thickness [m]
+    fields["FM"] = np.zeros((nrows, ncols))  # flow mass [kg]
+    fields["Vx"] = np.zeros((nrows, ncols))  # velocity in x direction [m/s]
+    fields["Vy"] = np.zeros((nrows, ncols))  # velocity in y direction [m/s]
+    fields["Vz"] = np.zeros((nrows, ncols))  # velocity in z direction [m/s]
+    fields["dmDet"] = np.zeros((nrows, ncols))  # flowing mass change due to detrainment [kg]
+    fields["FTStop"] = np.zeros((nrows, ncols))  # flow thickness that is stopped [m]
+    fields["FTDet"] = np.zeros((nrows, ncols))  # flow thickness that is detrained [m]
+    fields["FTEnt"] = np.zeros((nrows, ncols))  # flow thickness that is entrained [m]
+    fields["sfcChange"] = np.zeros((nrows, ncols))  # depth that changes the surface [m]
+    fields["sfcChangeTotal"] = np.zeros((nrows, ncols))  # total depth that changed the surface [m]
+    fields["demAdapted"] = np.zeros((nrows, ncols))  # adapted topography [m]
     # for optional fields, initialize with dummys (minimum size array). The cython functions then need something
     # even if it is empty to run properly
     if ("TA" in resTypesLast) or ("pta" in resTypesLast):
-        fields["pta"] = np.zeros((nrows, ncols))
-        fields["TA"] = np.zeros((nrows, ncols))
+        fields["pta"] = np.zeros((nrows, ncols))  # peak travel angle [°]
+        fields["TA"] = np.zeros((nrows, ncols))  # travel angle [°]
         fields["computeTA"] = True
         log.debug("Computing Travel Angle")
     else:
@@ -1681,15 +1693,15 @@ def initializeFields(cfg, dem, particles, releaseLine):
         fields["TA"] = np.zeros((1, 1))
         fields["computeTA"] = False
     if "pke" in resTypesLast:
-        fields["pke"] = np.zeros((nrows, ncols))
+        fields["pke"] = np.zeros((nrows, ncols))  # peak kinetic energy [kJ/m²]
         fields["computeKE"] = True
         log.debug("Computing Kinetic energy")
     else:
         fields["pke"] = np.zeros((1, 1))
         fields["computeKE"] = False
     if ("P" in resTypesLast) or ("ppr" in resTypesLast):
-        fields["P"] = np.zeros((nrows, ncols))
-        fields["ppr"] = np.zeros((nrows, ncols))
+        fields["P"] = np.zeros((nrows, ncols))  # pressure [kPa]
+        fields["ppr"] = np.zeros((nrows, ncols))  # peak pressure [kPa]
         fields["computeP"] = True
         log.debug("Computing Pressure")
     else:
@@ -2131,6 +2143,11 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     while t <= tEnd * (1.0 + 1.0e-13) and particles["iterate"]:
         startTime = time.time()
         log.debug("Computing time step t = %f s, dt = %f s" % (t, dt))
+
+        if cfgGen.getboolean("hydrograph"):
+            particles, fields, zPartArray0 = debF.releaseHydrograph(
+                cfg, inputSimLines, particles, fields, dem, zPartArray0, t
+            )
         # Perform computations
         particles, fields, zPartArray0, tCPU, dem = computeEulerTimeStep(
             cfgGen,
@@ -3594,6 +3611,7 @@ def adaptDEM(dem, fields, cfg):
     """adapt topography in respect to erosion and deposition
 
     Parameters
+    ---------
     dem: dict
         dictionary with info on DEM data
     fields : dict
diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini
index e00963570..121e0ef4d 100644
--- a/avaframe/com1DFA/com1DFACfg.ini
+++ b/avaframe/com1DFA/com1DFACfg.ini
@@ -132,6 +132,21 @@ entThDistVariation =
 # entrainment thickness (only considered if ENT file is shapefile and entThFromFile=False) [m]
 entTh =
 
+#+++++++++++++Hydrograph
+# if hydrograph is True, add hydrograph, provide the hydrograph values in a csv-table in the HYDR folder
+hydrograph = False
+# when checking if an already existing particle is within a hydrograph polygon
+# (that would cause numerical instabilities and rise an error), one can decide to add a buffer
+# around the polygon (0 means take strictly the points inside, a very small value
+# will include the points located on the polygon line)
+thresholdPointInHydr = 0.01
+# disabled at the moment:
+# distance that particles need to travel before new particles are initialized through the hydrograph
+# the distance is computed out of the hydrograph timesteps and the hydrograph velocity,
+#  if the velocity is 0, it's set to 1 m/s:
+# (distance = (timestep[i] - timestep[i-1]) * velocity)
+timeStepDistance = 5
+
 #++++++++++++Time stepping parameters
 # fixed time step (also used as first time step when using CFL) [s]
 dt = 0.1
diff --git a/avaframe/com1DFA/debrisFunctions.py b/avaframe/com1DFA/debrisFunctions.py
new file mode 100644
index 000000000..b0cb6eb96
--- /dev/null
+++ b/avaframe/com1DFA/debrisFunctions.py
@@ -0,0 +1,261 @@
+"""
+Functions that are used specifically for modeling debris flows (within DebrisFrame).
+"""
+
+# Load modules
+import logging
+import numpy as np
+import copy
+
+import avaframe.com1DFA.DFAfunctionsCython as DFAfunC
+import avaframe.com1DFA.particleTools as particleTools
+import avaframe.in3Utils.geoTrans as geoTrans
+import avaframe.com1DFA.com1DFA as com1DFA
+import avaframe.in2Trans.shpConversion as shpConv
+from avaframe.in1Data import getInput as gI
+
+# create local logger
+# change log level in calling module to DEBUG to see log messages
+log = logging.getLogger(__name__)
+
+
+def releaseHydrograph(cfg, inputSimLines, particles, fields, dem, zPartArray0, t, atol=1e-05):
+    """
+    Update particles with "new" particles initialised by a hydrograph.
+
+    Parameters
+    ---------
+    cfg: configparser
+        configuration settings
+    inputSimLines : dict
+        dictionary with input data dictionaries (releaseLine, hydrographAreaLine,...)
+    particles : dict
+        particles dictionary at t that are in the flow already
+    fields: dict
+        fields dictionary at t
+    dem: dict
+        dictionary with info on DEM data
+    zPartArray0: dict
+        dictionary containing z - value of particles at timestep 0
+    t: float
+        timestep of iteration
+    atol: float
+        look for matching time steps with atol tolerance - default is atol=1.e-5
+
+    Returns
+    ---------
+    particles: dict
+        particles dictionary at t including the hydrograph particles
+    fields: dict
+        updated fields dictionary at t including the hydrograph particles
+    zPartArray0: dict
+        dictionary containing z - value of particles at timestep 0
+    """
+    hydrValues = inputSimLines["hydrographAreaLine"]["values"]
+    if np.isclose(t, hydrValues["timeStep"], atol=atol, rtol=0).any():
+        iTup = np.where(np.isclose(t, hydrValues["timeStep"], atol=atol, rtol=0))
+        # iTup is a tuple containing an array with one value in the first position, so we can extract the index:
+        i = iTup[0].item()
+        log.info(
+            "add hydrograph at timestep: %f s with thickness %s m and velocity %s m/s"
+            % (t, hydrValues["thickness"][i], hydrValues["velocity"][i])
+        )
+        # similar workflow to secondary release!
+        particles, zPartArray0 = addHydrographParticles(
+            cfg,
+            particles,
+            inputSimLines,
+            hydrValues["thickness"][i],
+            hydrValues["velocity"][i],
+            dem,
+            zPartArray0,
+        )
+        particles = DFAfunC.getNeighborsC(particles, dem)
+        # update fields (compute grid values)
+        if fields["computeTA"]:
+            particles = DFAfunC.computeTrajectoryAngleC(particles, zPartArray0)
+        particles, fields = DFAfunC.updateFieldsC(cfg["GENERAL"], particles, dem, fields)
+
+    return particles, fields, zPartArray0
+
+
+def addHydrographParticles(cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0):
+    """
+    add new particles initialized by a hydrograph to particles that are in the flow already
+
+    Parameters
+    ---------
+    cfg: configparser object
+        configuration settings
+    particles : dict
+        particles dictionary at t that are in the flow already
+    inputSimLines : dict
+        dictionary with input data dictionaries (releaseLine, hydrographAreaLine,...)
+    thickness: float
+        thickness of incoming hydrograph
+    velocityMag: float
+        velocity of incoming hydrograph
+    dem: dict
+        dictionary with info on DEM data
+    zPartArray0: numpy array
+        z - value of particles at timestep 0
+
+    Returns
+    ---------
+    particles: dict
+        particles dictionary at t including the hydrograph particles
+    zPartArray0: dict
+        dictionary containing z - value of particles at timestep 0
+    """
+    hydrLine = inputSimLines["hydrographAreaLine"]
+    hydrLine["header"] = dem["originalHeader"].copy()
+    hydrLine = geoTrans.prepareArea(
+        hydrLine,
+        dem,
+        np.sqrt(2),
+        thList=[thickness],
+        combine=True,
+        checkOverlap=False,
+    )
+
+    # check if already existing particles are within the hydrograph polygon
+    # it's possible that there are still a few particles in the polygon with low velocities
+    # TODO: could think of a threshold of number of particles that are still allowed in the polygons?
+    mask = geoTrans.getParticlesInPolygon(
+        particles, hydrLine, cfg["GENERAL"].getfloat("thresholdPointInHydr")
+    )
+    if np.sum(mask) > 0:
+        # if there is at least one particle within the polygon (including the buffer):
+        message = (
+            "Already existing particles are within the hydrograph polygon, which can cause numerical instabilities (at timestep: %02f s)"
+            % (particles["t"] + particles["dt"])
+        )
+        # timestep in particles is not updated yet
+        log.error(message)
+        raise ValueError(message)
+
+    particlesHydrograph = com1DFA.initializeParticles(
+        cfg["GENERAL"],
+        hydrLine,
+        dem,
+    )
+    particlesHydrograph = DFAfunC.updateInitialVelocity(
+        cfg["GENERAL"], particlesHydrograph, dem, velocityMag
+    )
+
+    particles = particleTools.mergeParticleDict(particles, particlesHydrograph)
+    # save initial z position for travel angle computation
+    zPartArray0 = np.append(zPartArray0, copy.deepcopy(particlesHydrograph["z"]))
+    return particles, zPartArray0
+
+
+def checkHydrograph(hydrographValues, hydrCsv):
+    """
+    check if hydrograph satisfied the following requirements:
+    - hydrograph-timesteps are unique
+    - provided release-thickness is larger than zero
+    - the hydrograph-timesteps are not too close (that the particle density becomes too high)
+
+    Parameters
+    -----------
+    hydrCsv: str
+        directory to csv table containing hydrograph values
+    hydrographValues: dict
+        contains hydrograph values: timestep, thickness, velocity
+    """
+    # check if timesteps are unique
+    timeStepUnique = np.unique(hydrographValues["timeStep"])
+    if timeStepUnique.ndim == 0:
+        if timeStepUnique != hydrographValues["timeStep"]:
+            message = "The provided hydrograph time steps in %s are not unique" % (hydrCsv)
+            log.error(message)
+            raise ValueError(message)
+    elif len(timeStepUnique) != len(hydrographValues["timeStep"]):
+        message = "The provided hydrograph timesteps in %s are not unique" % (hydrCsv)
+        log.error(message)
+        raise ValueError(message)
+
+    # check that hydrograph thickness > 0
+    for th in hydrographValues["thickness"]:
+        if th <= 0:
+            message = "For every release time step a thickness > 0 needs to be provided in %s" % (hydrCsv)
+            log.error(message)
+            raise ValueError(message)
+
+
+def preparehydrographAreaLine(inputSimFiles, demOri, cfg):
+    """
+    read hydrograph polygon and values
+
+    Parameters
+    ----------
+    inputSimFiles : dict
+        dictionary containing
+        - hydrographFile: str, path to hydrograph polygon file
+        - hydrographCsv: str, path to hydrograph values (csv-)file
+    cfg: configparser object
+        configuration for simType
+    demOri : dict
+        dictionary with dem info (header original origin), raster data correct mesh cell size
+
+    Returns
+    -------
+    hydrLine: dict
+        contains hydrograph outline and values, among other things:
+        - x, y, z
+        - values: timeStep, thickness, velocity
+    """
+    try:
+        hydrFile = inputSimFiles["hydrographFile"]
+        hydrLine = shpConv.readLine(hydrFile, "", demOri)
+        hydrLine["fileName"] = hydrFile
+        hydrLine["type"] = "Hydrograph"
+        gI.checkForMultiplePartsShpArea(
+            cfg["GENERAL"]["avalancheDir"], hydrLine, "com1DFA", type="hydrograph"
+        )
+    except:
+        message = "No hydrograph shp file found"
+        log.error(message)
+        raise FileNotFoundError(message)
+
+    try:
+        hydrLine["values"] = gI.getHydrographCsv(inputSimFiles["hydrographCsv"])
+        hydrLine["thicknessSource"] = ["csv file"]
+    except:
+        message = "No hydrograph csv file found"
+        log.error(message)
+        raise FileNotFoundError(message)
+
+    checkHydrograph(hydrLine["values"], inputSimFiles["hydrographCsv"])
+
+    return hydrLine
+
+
+def checkTravelledDistance(cfgGen, hydrographValues, hydrCsv):
+    """
+    not used at the moment (related to timeStepDistance in the configuration file)!
+    check if time steps of hydrograph are not to close that the particle density becomes too high
+    check that particles moved out of hydrograph area before new particles are initialized
+    time between hydrograph time steps
+    first timestep is skipped since this is always ok.
+
+    Parameters
+    -----------
+    hydrCsv: str
+        directory to csv table containing hydrograph values
+    cfgGen: configparser
+        configuration settings, section "GENERAL"
+    hydrographValues: dict
+        contains hydrograph values: timestep, thickness, velocity
+    """
+    timeStepUnique = np.unique(hydrographValues["timeStep"])
+    if timeStepUnique.ndim > 0:
+        hydrDT = np.append(hydrographValues["timeStep"], 0) - np.append(0, hydrographValues["timeStep"])
+        vel = np.where(np.array(hydrographValues["velocity"]) > 0, np.array(hydrographValues["velocity"]), 1)
+        distance = vel[:-1] * hydrDT[1:-1]
+
+        if np.any(distance < cfgGen.getfloat("timeStepDistance")):
+            message = "Please select timesteps with greater spacing in %s." % (hydrCsv)
+            # TODO: error or warning?
+            log.error(message)
+            raise ValueError(message)
diff --git a/avaframe/in1Data/getInput.py b/avaframe/in1Data/getInput.py
index 1d165c9d7..40f7dc8db 100644
--- a/avaframe/in1Data/getInput.py
+++ b/avaframe/in1Data/getInput.py
@@ -296,6 +296,13 @@ def getInputDataCom1DFA(avaDir):
     entResInfo["resRemeshed"] = "No"
     entResInfo["bhdRemeshed"] = "No"
 
+    hydrographFile, entResInfo["hydrograph"], _ = getAndCheckInputFiles(
+        inputDir, "HYDR", "Hydrograph", fileExt="shp"
+    )
+    hydrographCsv, entResInfo["hydrographCsv"], _ = getAndCheckInputFiles(
+        inputDir, "HYDR", "Hydrograph parameters (csv)", fileExt="csv"
+    )
+
     # return DEM, first item of release, entrainment and resistance areas
     inputSimFiles = {
         "demFile": demFile,
@@ -310,6 +317,8 @@ def getInputDataCom1DFA(avaDir):
         "kFile": kFile,
         "tauCFile": tauCFile,
         "bhdFile": bhdFile,
+        "hydrographFile": hydrographFile,
+        "hydrographCsv": hydrographCsv,
     }
 
     for thFile in ["rel", "secondaryRel", "ent"]:
@@ -356,7 +365,7 @@ def getAndCheckInputFiles(inputDir, folder, inputType, fileExt="shp", fileSuffix
     """
     available = "No"
 
-    supportedFileFormats = [".shp", ".asc", ".tif"]
+    supportedFileFormats = [".shp", ".asc", ".tif", ".csv"]
 
     # Define the directory to search and the extensions
     if fileExt == "":
@@ -395,7 +404,8 @@ def getAndCheckInputFiles(inputDir, folder, inputType, fileExt="shp", fileSuffix
 
         if OutputFile.suffix not in supportedFileFormats:
             message = (
-                "Unsupported file format found for OutputFile %s; shp, asc, tif are allowed" % OutputFile
+                "Unsupported file format found for OutputFile %s; shp, asc, tif, csv are allowed"
+                % OutputFile
             )
             log.error(message)
             raise AssertionError(message)
@@ -1131,3 +1141,39 @@ def deriveLineRaster(
         log.info("%s read from %s" % (rasterNameStr[rasterType], str(rasterPath)))
 
     return rasterPath, lineDict
+
+
+def getHydrographCsv(hydrCsv):
+    """
+     get hydrograph values from the csv table
+     TODO: now the first column is defined as timestep, the second as thickness, third as velocity
+     see DebrisFrame Issue #18
+
+    Parameters
+    -----------
+    hydrCsv: str
+        path to csv file containing hydrograph values
+
+    Returns
+    -----------
+    hydrographValues: dict
+        contains hydrograph values: timestep, thickness, velocity
+    """
+    hydrParameters = np.genfromtxt(hydrCsv, delimiter=",", filling_values=0, skip_header=1)
+
+    if hydrParameters.ndim == 2:
+        # sort the columns according to the first column (timestep)
+        hydrParameters = hydrParameters[np.argsort(hydrParameters[:, 0])]
+        hydrographValues = {
+            "timeStep": hydrParameters[:, 0],
+            "thickness": hydrParameters[:, 1],
+            "velocity": hydrParameters[:, 2],
+        }
+    else:
+        hydrographValues = {
+            "timeStep": [hydrParameters[0]],
+            "thickness": [hydrParameters[1]],
+            "velocity": [hydrParameters[2]],
+        }
+
+    return hydrographValues
diff --git a/avaframe/in3Utils/geoTrans.py b/avaframe/in3Utils/geoTrans.py
index a2464d084..81f221dfe 100644
--- a/avaframe/in3Utils/geoTrans.py
+++ b/avaframe/in3Utils/geoTrans.py
@@ -1323,6 +1323,37 @@ def checkParticlesInRelease(particles, line, radius):
     particles : dict
         particles dictionary where particles outside of the polygon have been removed
     """
+    Mask = getParticlesInPolygon(particles, line, radius)
+    # also remove particles with negative mass
+    mask = np.where(particles["m"] <= 0, False, True)
+    Mask = np.logical_and(Mask, mask)
+    nRemove = len(Mask) - np.sum(Mask)
+    if nRemove > 0:
+        particles = particleTools.removePart(particles, Mask, nRemove, "")
+        log.debug("removed %s particles because they are not within the release polygon" % (nRemove))
+
+    return particles
+
+
+def getParticlesInPolygon(particles, line, radius):
+    """
+    check which particles are within a polygon (including a buffer)
+
+    Parameters
+    ----------
+    particles : dict
+        particles dictionary
+    line: dict
+        line dictionary of polygon
+    radius: float
+        this value is added to the polygon coordinates to decide whether a particles
+        is located within or outside of the polygon
+
+    Returns
+    -------
+    Mask: numpy array
+        is True at particle indices that are inside the polygon, and False for outside the polygon
+    """
     NameRel = line["Name"]
     StartRel = line["Start"]
     LengthRel = line["Length"]
@@ -1336,18 +1367,11 @@ def checkParticlesInRelease(particles, line, radius):
             "y": line["y"][int(start) : int(end)],
             "Name": name,
         }
+        # get an array that have at the same indices of a particle a True if the particle is within the polygon
+        # and a False if it is outside the polygon (including a buffer (radius))
         mask = pointInPolygon(line["header"], particles, avapath, radius)
         Mask = np.logical_or(Mask, mask)
-
-    # also remove particles with negative mass
-    mask = np.where(particles["m"] <= 0, False, True)
-    Mask = np.logical_and(Mask, mask)
-    nRemove = len(Mask) - np.sum(Mask)
-    if nRemove > 0:
-        particles = particleTools.removePart(particles, Mask, nRemove, "")
-        log.debug("removed %s particles because they are not within the release polygon" % (nRemove))
-
-    return particles
+    return Mask
 
 
 def pointInPolygon(demHeader, points, Line, radius):
@@ -1368,7 +1392,7 @@ def pointInPolygon(demHeader, points, Line, radius):
     Returns
     -------
     Mask : 1D numpy array
-        Mask of particles to keep
+        Mask of particles that are within the polygon
     """
     xllc = demHeader["xllcenter"]
     yllc = demHeader["yllcenter"]
diff --git a/avaframe/runStandardTestsCom1DFA.py b/avaframe/runStandardTestsCom1DFA.py
index 117020395..5c29c7c1e 100644
--- a/avaframe/runStandardTestsCom1DFA.py
+++ b/avaframe/runStandardTestsCom1DFA.py
@@ -254,10 +254,10 @@ def runSingleTest(
 # filter benchmarks for tag standardTest
 # filterType = "TAGS"
 # valuesList = ["resistance"]
-# filterType = "TAGS"
-# valuesList = ["standardTest", "standardTestSnowGlide"]
-filterType = "NAME"
-valuesList = ["avaInclinedPlaneEntresTest"]
+filterType = "TAGS"
+valuesList = ["standardTest", "standardTestSnowGlide"]
+# filterType = "NAME"
+# valuesList = ["avaInclinedPlaneEntresTest"]
 
 testList = tU.filterBenchmarks(testDictList, filterType, valuesList, condition="or")
 
diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py
index 488ccf91b..387f909b0 100644
--- a/avaframe/tests/test_DFAfunctionsCython.py
+++ b/avaframe/tests/test_DFAfunctionsCython.py
@@ -10,6 +10,7 @@
 import avaframe.com1DFA.DFAfunctionsCython as DFAfunC
 import avaframe.com1DFA.DFAtools as DFAtls
 import avaframe.com1DFA.DFAToolsCython as DFAtlsC
+import avaframe.in3Utils.geoTrans as geoTrans
 
 
 def test_getNeighborsC(capfd):
@@ -1005,6 +1006,72 @@ def test_computeForceC():
 
     assert (forceFrictArray == test_forceFrictArray).all()
 
+def test_updateInitialVelocity():
+    particles = {
+        "x": np.array([10.0]),
+        "y": np.array([10.0]),
+        "ux": np.array([0.0]),
+        "uy": np.array([0.0]),
+        "uz": np.array([0.0]),
+        "nPart": 1,
+        "velocityMag": np.array([0.0]),
+    }
+
+    # incline plane
+    dem = {
+        "header": {"nrows": 5, "ncols": 5, "cellsize": 5},
+        "rasterData": np.array(
+            [
+                [100, 100, 100, 100, 100],
+                [95, 95, 95, 95, 95],
+                [90, 90, 90, 90, 90],
+                [85, 85, 85, 85, 85],
+                [80, 80, 80, 80, 80],
+            ]
+        ),
+    }
+
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {"interpOption": "2"}
+
+    velocityMag = 10.0
+    dem = geoTrans.getNormalMesh(dem, num=1)
+
+    particlesTest = {
+        "x": np.array([10.0]),
+        "y": np.array([10.0]),
+        "ux": np.array([0.0]),
+        "uy": np.array([np.sqrt(50)]),
+        "uz": np.array([-np.sqrt(50)]),
+        "nPart": 1,
+        "velocityMag": np.array([velocityMag]),
+    }
+
+    particlesVelocity = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particles, dem, velocityMag)
+    for key in particlesTest:
+        assert np.isclose(particlesTest[key], particlesVelocity[key], atol=1e-4)
+    assert (
+        DFAtls.norm(particlesVelocity["ux"], particlesVelocity["uy"], particlesVelocity["uz"]) == velocityMag
+    )
+
+    # flat plane
+    dem["rasterData"] = np.ones((dem["header"]["ncols"], dem["header"]["nrows"]))
+    dem = geoTrans.getNormalMesh(dem, num=1)
+    particlesTest = {
+        "x": np.array([10.0]),
+        "y": np.array([10.0]),
+        "ux": np.array([0.0]),
+        "uy": np.array([0.0]),
+        "uz": np.array([0.0]),
+        "nPart": 1,
+        "velocityMag": np.array([0.0]),
+    }
+
+    particlesVelocity = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particles, dem, velocityMag)
+    for key in particlesTest:
+        assert np.isclose(particlesTest[key], particlesVelocity[key], atol=1e-3)
+
+
 """
 TODO: When calling pytest, the following function raises an error ("Fatal Python error: Aborted")
 (see issue #1002?)
diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py
index 6ac6f78cc..feec8bffd 100644
--- a/avaframe/tests/test_com1DFA.py
+++ b/avaframe/tests/test_com1DFA.py
@@ -1286,6 +1286,7 @@ def test_releaseSecRelArea():
         "thickness": [0.5, 1.0, 0.5],
         "rasterData": [secRelRaster1, secRelRaster2, secRelRaster3],
         "initializedFrom": "shapefile",
+        "type": "secondary release",
     }
     secondaryReleaseInfo["header"] = demHeader
     secondaryReleaseInfo["header"]["xllcenter"] = dem["originalHeader"]["xllcenter"]
@@ -1364,7 +1365,7 @@ def test_getRelThFromPart():
     # setup required input
     cfg = configparser.ConfigParser()
     cfg["GENERAL"] = {"relThFromFile": "True", "relTh": ""}
-    inputSimLines = {"releaseLine": {"thickness": ["1.2", "1.5"], "id": ["0", "1"]}}
+    inputSimLines = {"releaseLine": {"thickness": ["1.2", "1.5"], "id": ["0", "1"], "type": "Release"}}
     relThField = ""
 
     # call function to be tested
@@ -1441,6 +1442,7 @@ def test_initializeParticles():
         "Name": [""],
         "thickness": [1.0],
         "rasterData": relRaster,
+        "type": "Release",
     }
 
     releaseLine["header"] = demHeader
diff --git a/avaframe/tests/test_debrisFunctions.py b/avaframe/tests/test_debrisFunctions.py
new file mode 100644
index 000000000..0d6369069
--- /dev/null
+++ b/avaframe/tests/test_debrisFunctions.py
@@ -0,0 +1,298 @@
+""" Tests for module debrisFunctions """
+
+import numpy as np
+import pytest
+import configparser
+
+import avaframe.com1DFA.debrisFunctions as debF
+
+
+def test_addHydrographParticles():
+    inputSimLines = {
+        "hydrographAreaLine": {
+            "Name": ["testHydr"],
+            "Start": np.asarray([0.0]),
+            "Length": np.asarray([5]),
+            "type": "Hydrograph",
+            "x": np.asarray(
+                [
+                    0,
+                    10.0,
+                    10.0,
+                    0.0,
+                    0.0,
+                ]
+            )
+            - 2.5,
+            "y": np.asarray([0.0, 0.0, 10.0, 10.0, 0.0]) - 2.5,
+            "thicknessSource": ["csv file"],
+            "thickness": 1,
+        }
+    }
+    thickness = inputSimLines["hydrographAreaLine"]["thickness"]
+    velocityMag = 0
+
+    demHeader = {}
+    demHeader["xllcenter"] = 0
+    demHeader["yllcenter"] = 0
+    demHeader["cellsize"] = 5.0
+    demHeader["nodata_value"] = -9999
+    demHeader["nrows"] = 7
+    demHeader["ncols"] = 7
+    dem = {"header": demHeader}
+    dem["rasterData"] = np.ones((demHeader["nrows"], demHeader["ncols"]))
+    dem["originalHeader"] = dem["header"]
+    dem["areaRaster"] = np.ones((demHeader["nrows"], demHeader["ncols"]))
+    dem["Nx"] = np.zeros_like(dem["rasterData"])
+    dem["Ny"] = np.zeros_like(dem["rasterData"])
+    dem["Nz"] = np.zeros_like(dem["rasterData"])
+
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {
+        "resType": "ppr|pft|pfv",
+        "rho": "1000.",
+        "gravAcc": "9.81",
+        "cpIce": "2050",
+        "TIni": "-10",
+        "avalancheDir": "data/avaKotHYDR",
+        "massPerParticleDeterminationMethod": "MPPDH",
+        "interpOption": "2",
+        "initialiseParticlesFromFile": "False",
+        "iniStep": "False",
+        "seed": "12345",
+        "sphKernelRadius": "1",
+        "deltaTh": "1",
+        "initPartDistType": "uniform",
+        "thresholdPointInPoly": "0.001",
+        "massPerPart": "1000",
+        "thresholdPointInHydr": "0",
+    }
+
+    particles = {
+        "nPart": 3,
+        "x": np.array([12, 20, 30]),
+        "y": np.array([5, 10, 30]),
+        "z": np.array([1, 1, 1]),
+        "m": np.array([1000, 1000, 1000]),
+        "idFixed": np.array([0, 0, 0]),
+        "t": 1.0,
+        "dt": 0.1,
+    }
+    nPart = particles["nPart"]
+    particles["totalEnthalpy"] = (
+        cfg["GENERAL"].getfloat("TIni") * cfg["GENERAL"].getfloat("cpIce")
+        + cfg["GENERAL"].getfloat("gravAcc") * particles["z"]
+    )
+    particles["massPerPart"] = 1000
+    particles["mTot"] = np.sum(particles["m"])
+    particles["tPlot"] = 0
+    particles["h"] = np.ones(nPart)
+    particles["ux"] = np.zeros(nPart)
+    particles["uy"] = np.zeros(nPart)
+    particles["uz"] = np.zeros(nPart)
+    particles["uAcc"] = np.zeros(nPart)
+    particles["velocityMag"] = np.zeros(nPart)
+    particles["trajectoryLengthXY"] = np.zeros(nPart)
+    particles["trajectoryLengthXYCor"] = np.zeros(nPart)
+    particles["trajectoryLengthXYZ"] = np.zeros(nPart)
+    particles["trajectoryAngle"] = np.zeros(nPart)
+    particles["stoppCriteria"] = False
+    particles["peakForceSPH"] = 0.0
+    particles["forceSPHIni"] = 0.0
+    particles["peakMassFlowing"] = 0
+    particles["xllcenter"] = dem["originalHeader"]["xllcenter"]
+    particles["yllcenter"] = dem["originalHeader"]["yllcenter"]
+    particles["nExitedParticles"] = 0.0
+    particles["dmDet"] = np.zeros(nPart)
+    particles["dmEnt"] = np.zeros(nPart)
+
+    zPartArray0 = np.array([1, 1, 1])
+    newParticleNumber = 4
+    particlesTest = {
+        "nPart": newParticleNumber + 3,
+        "mTot": 7000,
+        "x": np.append(particles["x"], np.array([0, 5, 0, 5])),
+        "y": np.append(particles["y"], np.array([0, 0, 5, 5])),
+        "z": np.append(particles["z"], np.ones([newParticleNumber])),
+        "m": np.append(particles["m"], np.ones([newParticleNumber]) * 1000),
+    }
+    zPartArray0Test = np.ones(particlesTest["nPart"])
+
+    particlesHydr, zPartArray0Hydr = debF.addHydrographParticles(
+        cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0
+    )
+
+    assert np.all(np.equal(zPartArray0Hydr, zPartArray0Test))
+    for key in particlesTest:
+        if key in ["nPart", "mTot"]:
+            assert particlesTest[key] == particlesHydr[key]
+        else:
+            assert np.all(np.equal(particlesTest[key], particlesHydr[key]))
+    for key in ["ux", "uy", "uz", "velocityMag"]:
+        assert np.all(np.equal(np.zeros(particlesTest["nPart"]), particlesHydr[key]))
+
+    cfg["GENERAL"]["deltaTh"] = "0.25"
+    cfg["GENERAL"]["initPartDistType"] = "random"
+    cfg["GENERAL"]["thresholdMassSplit"] = "1.5"
+
+    particlesHydr, zPartArray0Hydr = debF.addHydrographParticles(
+        cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0
+    )
+    assert particlesHydr["nPart"] == 16 + 3
+    for key in ["ux", "uy", "uz", "velocityMag", "x", "y", "z"]:
+        assert len(particlesHydr[key]) == particlesHydr["nPart"]
+    assert particlesHydr["mTot"] == 7000
+
+    particles["x"] = np.array([4, 10, 30])
+    particles["y"] = np.array([5, 3, 30])
+
+    with pytest.raises(ValueError):
+        debF.addHydrographParticles(cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0)
+
+"""
+TODO: When calling pytest, the following function raises an error ("Fatal Python error: Aborted")
+(see issue #1002?)
+
+def test_releaseHydrograph():
+    inputSimLines = {
+        "hydrographAreaLine": {
+            "Name": ["testHydr"],
+            "Start": np.asarray([0.0]),
+            "Length": np.asarray([5]),
+            "type": "Hydrograph",
+            "x": np.asarray(
+                [
+                    0,
+                    10.0,
+                    10.0,
+                    0.0,
+                    0.0,
+                ]
+            )
+            - 2.5,
+            "y": np.asarray([0.0, 0.0, 10.0, 10.0, 0.0]) - 2.5,
+            "thicknessSource": ["csv file"],
+            "thickness": 1,
+            "values": {
+                "timeStep": np.array([0, 1]),
+                "thickness": np.array([1, 1]),
+                "velocity": np.array([0, 0]),
+            },
+        }
+    }
+    thickness = inputSimLines["hydrographAreaLine"]["thickness"]
+    velocityMag = 0
+
+    demHeader = {}
+    demHeader["xllcenter"] = 0
+    demHeader["yllcenter"] = 0
+    demHeader["cellsize"] = 5.0
+    demHeader["nodata_value"] = -9999
+    demHeader["nrows"] = 7
+    demHeader["ncols"] = 7
+    dem = {"header": demHeader}
+    dem["rasterData"] = np.ones((demHeader["nrows"], demHeader["ncols"]))
+    dem["originalHeader"] = dem["header"]
+    dem["areaRaster"] = np.ones((demHeader["nrows"], demHeader["ncols"]))
+    dem["Nx"] = np.zeros_like(dem["rasterData"])
+    dem["Ny"] = np.zeros_like(dem["rasterData"])
+    dem["Nz"] = np.zeros_like(dem["rasterData"])
+    dem["headerNeighbourGrid"] = demHeader
+
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {
+        "resType": "ppr|pft|pfv",
+        "rho": "1000.",
+        "rhoEnt": "1000.",
+        "gravAcc": "9.81",
+        "cpIce": "2050",
+        "TIni": "-10",
+        "avalancheDir": "data/avaKotHYDR",
+        "massPerParticleDeterminationMethod": "MPPDH",
+        "interpOption": "2",
+        "initialiseParticlesFromFile": "False",
+        "iniStep": "False",
+        "seed": "12345",
+        "sphKernelRadius": "1",
+        "deltaTh": "1",
+        "initPartDistType": "uniform",
+        "thresholdPointInPoly": "0.001",
+        "massPerPart": "1000",
+    }
+
+    particles = {
+        "nPart": 3,
+        "x": np.array([10, 20, 30]),
+        "y": np.array([5, 10, 30]),
+        "z": np.array([1, 1, 1]),
+        "m": np.array([1000, 1000, 1000]),
+        "idFixed": np.array([0, 0, 0]),
+        "dmDet": np.array([0, 0, 0]),
+        "dmEnt": np.array([0, 0, 0]),
+        "ux": np.array([0, 0, 0]),
+        "uy": np.array([0, 0, 0]),
+        "uz": np.array([0, 0, 0]),
+        "trajectoryAngle": np.array([0, 0, 0]),
+        "stoppedParticles": {
+            "m": np.empty(0),
+            "x": np.empty(0),
+            "y": np.empty(0),
+        },
+    }
+    nPart = particles["nPart"]
+    particles["totalEnthalpy"] = (
+        cfg["GENERAL"].getfloat("TIni") * cfg["GENERAL"].getfloat("cpIce")
+        + cfg["GENERAL"].getfloat("gravAcc") * particles["z"]
+    )
+    particles["massPerPart"] = 1000
+    particles["mTot"] = np.sum(particles["m"])
+    particles["tPlot"] = 0
+    particles["h"] = np.ones(nPart)
+    particles["uAcc"] = np.zeros(nPart)
+    particles["velocityMag"] = np.zeros(nPart)
+    particles["trajectoryLengthXY"] = np.zeros(nPart)
+    particles["trajectoryLengthXYCor"] = np.zeros(nPart)
+    particles["trajectoryLengthXYZ"] = np.zeros(nPart)
+    particles["stoppCriteria"] = False
+    particles["peakForceSPH"] = 0.0
+    particles["forceSPHIni"] = 0.0
+    particles["peakMassFlowing"] = 0
+    particles["xllcenter"] = dem["originalHeader"]["xllcenter"]
+    particles["yllcenter"] = dem["originalHeader"]["yllcenter"]
+    particles["nExitedParticles"] = 0.0
+
+    fields = {
+        "computeTA": 0,
+        "computeKE": 0,
+        "computeP": 0,
+        "pfv": np.zeros_like(dem["rasterData"]),
+        "ppr": np.zeros_like(dem["rasterData"]),
+        "pft": np.zeros_like(dem["rasterData"]),
+        "pta": np.zeros_like(dem["rasterData"]),
+        "pke": np.zeros_like(dem["rasterData"]),
+        "dmDet": np.zeros_like(dem["rasterData"]),
+    }
+    fields["pft"][[2, 4, 6], [1, 2, 6]] = 1
+
+    zPartArray0 = np.array([1, 1, 1])
+
+    t = 0.00
+
+    newParticleNumber = 4
+    particlesTest = {
+        "nPart": newParticleNumber + 3,
+        "mTot": 7000,
+        "x": np.append(particles["x"], np.array([0, 5, 0, 5])),
+        "y": np.append(particles["y"], np.array([0, 0, 5, 5])),
+        "z": np.append(particles["z"], np.ones([newParticleNumber])),
+        "m": np.append(particles["m"], np.ones([newParticleNumber]) * 1000),
+    }
+    zPartArray0Test = np.ones(particlesTest["nPart"])
+
+    debF.updateParticlesHydrograph(cfg, inputSimLines, particles, fields, dem, zPartArray0, t)
+
+"""
+
+
+if __name__ == "__main__":
+    test_addHydrographParticles()
diff --git a/avaframe/tests/test_geoTrans.py b/avaframe/tests/test_geoTrans.py
index 30bdad6ba..6f20facad 100644
--- a/avaframe/tests/test_geoTrans.py
+++ b/avaframe/tests/test_geoTrans.py
@@ -510,6 +510,43 @@ def test_checkParticlesInRelease():
     assert particles2["nPart"] == 4
 
 
+def test_getParticlesInPolygon():
+    """test if particles are within release polygon"""
+    # setup required input
+    hydrographLine = {
+        "Name": ["testRel"],
+        "Start": np.asarray([0.0]),
+        "Length": np.asarray([5]),
+        "type": "Release",
+        "x": np.asarray([0, 10.0, 10.0, 0.0, 0.0]),
+        "y": np.asarray([0.0, 0.0, 10.0, 10.0, 0.0]),
+    }
+    demHeader = {}
+    demHeader["xllcenter"] = 0.0
+    demHeader["yllcenter"] = 0.0
+    demHeader["cellsize"] = 1.0
+    demHeader["noDataValue"] = -9999
+    demHeader["nrows"] = 5
+    demHeader["ncols"] = 5
+    hydrographLine["header"] = demHeader
+    particles = {
+        "x": np.asarray([2.4, 9.7, 10.02, 11.5]),
+        "y": np.asarray([2.4, 9.7, 10.2, 11.5]),
+        "nPart": 4,
+        "m": np.asarray([1.4, 1.7, 1.4, 1.8]),
+    }
+    radius = 0
+
+    MaskTest = np.array([True, True, False, False])
+    Mask = geoTrans.getParticlesInPolygon(particles, hydrographLine, radius)
+    assert np.array_equal(MaskTest, Mask)
+
+    radius = 0.5
+    MaskTest = np.array([True, True, True, False])
+    Mask = geoTrans.getParticlesInPolygon(particles, hydrographLine, radius)
+    assert np.array_equal(MaskTest, Mask)
+
+
 def test_remeshDEM(tmp_path):
     """test size of interpolated data onto new mesh"""
 
diff --git a/avaframe/tests/test_particleInitialisation.py b/avaframe/tests/test_particleInitialisation.py
index dbb631afa..bd5c51b64 100644
--- a/avaframe/tests/test_particleInitialisation.py
+++ b/avaframe/tests/test_particleInitialisation.py
@@ -163,6 +163,7 @@ def test_getIniPosition(tmp_path):
             "header": demOri["header"],
             "y": np.asarray([5.0, 5.0, 10.0, 10.0, 5.0]),
             "rasterData": relRaster,
+            "type": "Release",
         }
     }