diff --git a/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini b/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini index 852e8fc12..423307056 100644 --- a/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini +++ b/avaframe/com6RockAvalanche/com6RockAvalancheCfg.ini @@ -48,3 +48,64 @@ frictModel = Voellmy #+++++++++++++Voellmy friction model muvoellmy = 0.035 xsivoellmy = 700. +# expected mesh size [m] +meshCellSize = 5 + +# density of entrained snow [kg/m³] +rhoEnt = 100 +#+++++Entrainment thickness++++ +# True if entrainment thickness should be read from file (shapefile, raster); if False - entTh read from ini file (only available for shapefile) +entThFromFile = True +# if a thickness value is missing for the entrainment feature in the provided shp file this value is used for all features [m] +entThIfMissingInShp = 0.3 +# VARIATION options only if ENT file is a shapefile -------- +# if a variation on entTh shall be performed add here +- percent and number of steps separated by $ +# for example entThPercentVariation=50$10 [%] +entThPercentVariation = +# if a variation on entTh shall be performed add here +- absolute value and number of steps separated by $ +# for example entThRangeVariation=0.5$10 [m] +entThRangeVariation = +# if a variation on entTh shall be performed add here +- ci% value and number of steps separated by $ +# for example entThRangeFromCiVariation= ci95$10 +entThRangeFromCiVariation = +# if variation on entTh shall be performed using a normal distribution in number of steps, +# value of buildType (ci95 value), min and max of dist in percent, buildType (ci95 only allowed), +# support (e.g. 10000) all separated by $: e.g. normaldistribution$numberOfSteps$0.3$95$ci95$10000 +# if entFromShp=True ci95 is read from shp file too +entThDistVariation = +# entrainment thickness (only considered if ENT file is shapefile and entThFromFile=False) [m] +entTh = + + +#++++++++++++ Entrainment Erosion Energy +# Used to determine speed loss via energy loss due to entrained mass +entEroEnergy = 5000 +entShearResistance = 0 +entDefResistance = 0 + +#++++++++++++ Resistance model +# default setup: +ResistanceModel = default +# At each time step, ResistanceModel default applies increased friction and optional detrainment (see below). Only relevant in resistance areas. +# NOTE: development setup; parameter values need more testing and calibration!! + +# parameter for increased friction in resistance areas +cResH = 0.01 + +# Apply detrainment in resistance areas in accordance with the flow thickness and flow velocity thresholds specified below. +# if False - only increased friction is applied in resistance areas +detrainment = True + +# detrainment parameter defined by Feistl et al. (2014) +detK = 5 + +# thresholds if detrainment is set to True +# FV OR FT below min thresholds: apply only detrainment. no increased friction +# FV AND FT within min and max thresholds: no detrainment, only apply increased friction +# FV OR FT above max thresholds: no detrainment and no increased friction applied +forestVMin = 6. +forestThMin = 0.6 +forestVMax = 40. +forestThMax = 10. + + diff --git a/avaframe/com6RockAvalanche/scarp.py b/avaframe/com6RockAvalanche/scarp.py index d3ae2f762..dc3b9b938 100644 --- a/avaframe/com6RockAvalanche/scarp.py +++ b/avaframe/com6RockAvalanche/scarp.py @@ -149,8 +149,14 @@ def scarpAnalysisMain(cfg, baseDir): ) else: raise ValueError("Unsupported method. Choose 'plane' or 'ellipsoid'.") + hRelData = dem["rasterData"] - scarpData + min_hrel = hRelData.min() + if min_hrel < 0: + log.warning("hRel: most negative clipped value = %.2f m", min_hrel) + hRelData = np.maximum(0, hRelData) #enforce non-negative release + #Compute and log excavated volume cellArea = abs(dem["header"]["cellsize"] ** 2) @@ -247,6 +253,8 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes): betaX = [ math.tan(slopeRad) * math.sin(dipRad) ] betaY = [ math.tan(slopeRad) * math.cos(dipRad) ] + min_clipped = 0.0 + for i in range(1, nPlanes): xSeed.append(planes[5 * i]) ySeed.append(planes[5 * i + 1]) @@ -266,11 +274,17 @@ def calculateScarpWithPlanes(elevData, periData, elevTransform, planes): scarpVal = zSeed[0] + (north - ySeed[0]) * betaY[0] - (west - xSeed[0]) * betaX[0] for k in range(1, nPlanes): scarpVal = max(scarpVal, zSeed[k] + (north - ySeed[k]) * betaY[k] - (west - xSeed[k]) * betaX[k]) - + if periData[row, col] > 0: - scarpData[row, col] = min(elevData[row, col], scarpVal) + val = min(elevData[row, col], scarpVal) + if val < 0: + min_clipped = min(min_clipped, val) + scarpData[row, col] = max(0, val) else: scarpData[row, col] = elevData[row, col] + + if min_clipped < 0: + log.warning("Plane scarp: most negative clipped value = %.2f m", min_clipped) return scarpData @@ -305,6 +319,7 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids): xCenter, yCenter, maxDepth = [], [], [] semiMajor, semiMinor = [], [] tilt, tiltDir, offset, dip = [], [], [], [] + min_clipped = 0.0 for i in range(nEllipsoids): xCenter.append(ellipsoids[9 * i]) @@ -373,6 +388,15 @@ def calculateScarpWithEllipsoids(elevData, periData, elevTransform, ellipsoids): totalDepth = baseDepth + tiltEffect + z0 scarpVal = min(scarpVal, elevData[row, col] - totalDepth) - scarpData[row, col] = scarpVal if periData[row, col] > 0 else elevData[row, col] + if periData[row, col] > 0: + val = min(elevData[row, col], scarpVal) + if val < 0: + min_clipped = min(min_clipped, val) + scarpData[row, col] = max(0, val) + else: + scarpData[row, col] = elevData[row, col] + + if min_clipped < 0: + log.warning("Ellipsoid scarp: most negative clipped value = %.2f m", min_clipped) return scarpData