Methodology

What the simulator computes, what it assumes, how each part was verified, and where it stops being valid.

Scope and intended use

This is a screening and design-support tool for rockfall hazard over a drone-derived DTM. It produces hazard rasters, runout statistics, barrier design values, and supporting terrain analyses.

Terrain: resolution, units and nodata

Resolution

Survey DTMs are routinely gigapixel — the reference dataset is 29 702 × 53 127 float32, 6.3 GB for the full-resolution band alone. Reading that band would need tens of gigabytes of RAM, so the loader instead reads from the file's overview pyramid, choosing the cheapest level that gets within a quarter of the requested output size, and decimating from there. A GeoTIFF with no pyramid is streamed strip by strip and decimated on the fly, so peak memory stays bounded either way.

The consequence is the single most important thing to understand about every number the app produces: the simulation runs on the decimated grid, not on the native survey. The cell size actually used is reported in the results table and written into every export.

Units

The physics is Newtonian and works in metres. A DTM in a geographic CRS has its pixel size expressed in degrees — for the reference dataset 1.077 × 10⁻⁶°, which taken literally is a one-micrometre cell and makes every velocity, energy and runout meaningless. Geographic grids are therefore converted on load to a local east/north metric frame using the exact metres-per-degree at the grid's centre latitude (accurate to millimetres over a site-sized extent).

A grid with square pixels in degrees is not square in metres. The reference DTM's 5.24 m × 5.73 m cells are carried through as two distinct dimensions rather than averaged; the eight D8 neighbour distances, the mesh geometry and the cell area all use the real values.

Nodata

Two conventions have to be handled. The declared value from the GDAL_NODATA tag is honoured directly. Beyond that, exports frequently write a constant into unsurveyed ground without declaring it, and left alone that renders as a flat plane and enters every statistic as real terrain. Such a fill value is identified as a spike in the elevation histogram that sits at an extreme of the range and is separated from real ground by a gap no surface would leave — on the second reference DTM, 0.0 across 90% of the grid with the nearest real elevation 336 m above it. Ground that genuinely reaches 0 m rises continuously out of it, leaves no gap, and is not flagged; the value can also be set explicitly.

Nodata is honoured everywhere rather than interpolated across. The reference survey covers only 38% of its bounding box, and the uncovered area is a hole in the 3D mesh, a boundary where flow leaves the hydrological domain, and ground on which no rock can be released.

Imagery

Both imagery sources are resampled onto the DTM's own grid rather than stretched to fit its bounding box, so the texture registers with the terrain cell for cell.

An orthophoto is mapped through its own geotransform into the DTM's source CRS; the two must match, which is checked. A web basemap needs the DTM located in longitude/latitude, so a projected survey is converted with an inverse transverse Mercator (Snyder). That implementation round-trips against an independent forward transform to better than a millimetre in both hemispheres, and the web-Mercator tile arithmetic reproduces known tile corner positions to 4 × 10⁻¹³ pixels — the two together are the whole georeferencing chain.

Zoom is chosen as the deepest level whose tile count stays within budget, so tile resolution is finer than the DTM cell size and the drape is not the limiting factor. The visual check that matters is continuity: linear features crossing the survey boundary onto the surrounding context plane should not step.

Hillshade over imagery uses an overlay composite rather than flat alpha. Flat alpha mixes grey into every pixel and desaturates; overlay leaves mid-grey neutral and only darkens or lightens, which is what shaded relief is meant to do.

The trajectory model

Each rock is a point mass integrated with a fixed time step. One loop produces two regimes:

• Flight — parabolic projectile motion while the rock is above the surface.
• Contact — velocity is split into components normal and tangential to the surface, and one of two things happens.

The distinction inside contact is the part that most often goes wrong in simple implementations:

• A genuine impact (the rock was airborne and has now met the ground) reflects the normal velocity scaled by R_n and scales the tangential velocity by R_t, once.
• Sustained contact is treated as rolling/sliding: the into-surface normal component is removed and the tangential motion decelerates under Coulomb friction only.

Gravity acts every step. A rock stops when it stays below stop_speed while in contact for stop_hold consecutive steps. Outputs are aggregated onto the DTM grid as passing frequency, peak velocity, and peak translational kinetic energy (½mv²) — the last is derived from the velocity raster, since energy is monotonic in speed.

Release model: samples or blocks

rocks_per_source means two different things, and conflating them is a real error:

Release points are scattered within their cell rather than placed at the centre; identical start points make rocks from one cell follow near-identical paths. In swarm mode the scatter additionally enforces a minimum separation so blocks do not start inside one another, falling back to valid ground in preference to ideal spacing near the survey boundary.

Rock–rock collisions

All rocks advance on one shared clock. Blocks are spheres of the equivalent-volume radius. Contacts are found by sweep-and-prune along the widest-spread axis and resolved with an impulse between equal masses, using a normal restitution and a tangential damping term.

Resolution iterates (projected Gauss-Seidel): separating one pair can drive a rock into a third, and rocks funnelling into the same gully form exactly such piles. A single sweep left blocks interpenetrating by up to 43% of a block diameter; eight passes reduce that to under 0.2%. The initial packing is relaxed the same way before the run starts, alternating separation with re-seating on the surface, because on a steep face those two fight each other.

Post-event terrain

Each release cell loses one block's volume; each rock adds its volume where it stops. The result is mass-conserving apart from blocks that run off the edge of the DTM, which are reported separately.

Scar depth is exactly block volume ÷ cell area. On a 5.5 m grid a 0.30 m³ block is a 1 cm change — correct, but invisible as geometry, which is why the displayed post-event surface has a clearly labelled exaggeration control. Elevation change is averaged over each mesh cell rather than point-sampled, so the volume moved survives mesh decimation even when a scar is only one or two cells wide.

Barriers

A barrier is used in two modes.

As a measuring line it does not stop anything: each rock's first crossing is recorded with its energy and its height above ground. Sizing a fence from an unobstructed run is the correct order of operations — the design question is what would arrive if nothing were there. Design values are the 95th percentile, with 0.5 m of freeboard on the height.

As an obstacle the fence acts during the run, with three outcomes: retained (below the crest, within capacity), overtopped (airborne above the crest), or breached (energy exceeds capacity; the rock continues with the rated capacity removed from its kinetic energy). That last is a first-order treatment — the structure absorbs what it is rated for and fails — not a structural model of the fence.

Change detection

Standard repeat-survey geomorphic change detection, in three stages.

Co-registration

A residual horizontal shift between two flights appears in the difference as a slope- and aspect-dependent pattern. A linearised fit against the terrain gradients — equivalent to the published cosine form but without dividing by tan(slope) — supplies a starting estimate, which is then refined by direct pattern search on the actual resampling operator. The linearisation systematically misjudges the shift on rough ground because it uses central-difference gradients; searching the real objective has no such bias.

The correction is split half onto each surface. Interpolation smooths, and smoothing only one side makes rough ground differ everywhere.

Level of detection

The threshold is fitted per cell as σ = a + b·tan(slope) + c·roughness. A single flat threshold flags every steep face as change: residual misregistration becomes an elevation difference of shift × slope, and resampling error scales with how much the surface curves inside one cell. The fit is on absolute deviations with outlier rejection, so genuine change does not inflate the threshold meant to detect it.

Segmentation

Volumes are integrated over connected features that clear both the threshold and a minimum area. Summing every cell that individually clears a 95% threshold is what makes a noisy pair report thousands of cubic metres of nothing — one cell in twenty crosses it by chance, and across a million cells that is a large fake volume.

Calibration and identifiability

The objective is spatial, not volumetric: for every deposit, observed or simulated, the distance to the nearest release cell, compared as a full distribution via a 1-D Wasserstein distance over matched quantiles. This follows directly from the finding above — fitting to a measured volume would fit its geometry-dependent bias. Using the whole distribution rather than a summary statistic also penalises a model that gets the mean right by pairing too many short runs with too many long ones.

Fitting three parameters to one distribution need not have a unique answer, so the search keeps the entire objective grid and reports, per parameter, the range that stays within the seed-to-seed noise floor of the best fit — the misfit between two independent realisations of the same parameters. Anything closer than that is indistinguishable.

Result on the reference site

The noise floor is 0.192 — larger than the entire variation in the R_n and R_t profiles. Runout distance constrains rolling friction and says essentially nothing about restitution. That is consistent with the independent finding that rocks are airborne only ~2% of the time here: restitution acts only at impacts during flight, so if blocks roll nearly the whole way, friction does all the work.

Practically: constrain μ from runout, take R_n/R_t from literature for the rock type, and note the results are insensitive to them on this terrain anyway. Constraining restitution needs observations that depend on it — bounce heights or impact positions — which needs finer terrain detail than a 5 m grid.

Hydrology

Sinks are filled by priority flood, processing cells in order of increasing elevation inward from the domain edge, which guarantees each cell ends at least as high as its lowest path out in a single pass. Flow routing is D8 with the eight neighbour distances computed from the real non-square cell dimensions. Accumulation walks the grid in descending filled elevation, which is a valid topological order for a D8 network by construction. A picked outlet is snapped to the largest nearby channel before the catchment is traced upstream.

After filling, pit floors are flat and have no downhill neighbour, so flow terminates there rather than continuing — 5% of cells on the reference DTM. Mass is still conserved (the water stops rather than vanishing) but streams can end early in filled flats.

Objects, exposure and obstruction

The nDSM (DSM minus DTM) is the height of everything standing on the ground, and is the one product neither surface gives alone. The pair is co-registered first with the same estimator the change detection uses, because a residual horizontal shift becomes an apparent object height of shift × slope on every slope in the scene.

Objects are separated from ground by height, then classified by the coefficient of variation of object height within each object. A building has a coherent top, so its height above ground is near-constant even where it sits on a slope; a canopy is broken, so its variation is a large fraction of its height. An earlier version used a local 3×3 roughness window and failed: for an object a few cells across, most cells are edge cells where height jumps from zero to full building height, so every small building looked as rough as a tree.

Obstruction gives the two classes different physics, because they do different things. A building is rigid: a block arriving below roof height stops, and the energy delivered is recorded — that is the number a risk assessment needs. Vegetation is compliant: it does not arrest a block, it bleeds speed exponentially with the distance travelled inside the canopy at a stated rate. Verified on a test belt across the runout: buildings arrested 49% of blocks and cut mean runout from 75.3 m to 35.4 m, while an equivalent vegetation belt arrested none and cut runout to 69.7 m.

Running it as a service

Two assumptions the local app was built on stop holding once it is shared, and both were corrected rather than worked around.

State was global. One dictionary held the loaded terrain and its results, which is right for one person on one machine. Shared, it meant a second user's load replaced the first user's terrain, and the first user's next request computed against the wrong surface — with no error raised. State is now keyed by session. This is also why the service runs as a single instance: the sessions live in its memory, and a second instance would recreate the same fault across processes.

The DTM is no longer local. Rasters live in a bucket in the same region as the compute, and are read by byte range rather than downloaded. That works because the reader already streams pyramid overviews strip by strip: loading the 8.42 GB reference DTM issues about 253 ranged requests. The decimated grid is then cached at 1.4 MB, so the expensive read happens once per file rather than once per process.

What a saved run contains

A saved run holds the parameters, the results and the identity of the DTM it was computed on. All three are needed: parameters alone cannot be checked against anything, results alone cannot be explained, and either without the terrain could be drawn over the wrong surface.

The terrain itself is not stored. It is large, it is already cached, and a copy per run would multiply gigabytes for nothing; the run records the file, the grid resolution and the nodata handling needed to rebuild it. On recall the DTM is reloaded if it differs, and a run whose rasters do not match the current grid is refused rather than drawn — a hazard raster is defined cell by cell against one particular surface.

Serialisation walks the result structure rather than naming fields. Analyses return nested dictionaries of arrays and those structures change as the tool grows; a hand-maintained field list would quietly stop capturing new results, and a run that silently lost half its contents is worse than one that failed to save. Arrays keep their exact dtype and their NaNs, and rendered images are stored byte for byte, so a recalled run is identical rather than equivalent.

Flood screening

Inundation is mapped with HAND — height above nearest drainage (Rennó et al. 2008; Nobre et al. 2011). Each cell is referenced to the channel cell it actually drains to along the D8 network, so its height is measured along the flow path rather than perpendicular to the nearest stream. Inundation at a channel stage is then HAND ≤ stage and depth is stage − HAND. This is the basis of operational first-pass inundation mapping.

HAND is computed in a single ascending-elevation pass: a cell's downstream neighbour is always lower, so by the time a cell is visited its downstream reference already exists and is simply inherited. Verified against a brute-force walk of every cell's flow path to its channel — identical to 0.00e+00 on three different surfaces, with identical undefined sets.

In rainfall mode, peak discharge comes from the rational method (Q = C i A, which assumes the whole contributing area responds at once — reasonable for small steep catchments, over-conservative for large ones) and stage from Manning's equation for a wide rectangular channel, d = (Qn / (W√S))3/5. Channel width is an input, not measured, because a 16 m grid cannot resolve a wadi cross-section. Verified: discharge back-substitutes through Manning to the input value, and depth scales as Q^0.6 exactly.

Slope kinematics

Markland tests ask a purely geometric question: given the face orientation, this joint set and a friction angle φ, is failure kinematically possible?

• Planar — the joint daylights in the face (dip directions within tolerance), is shallower than the face, and is steeper than friction: face dip > joint dip > φ.
• Wedge — the intersection line of two joints, from the cross product of their normals, plunges less than the face and more than φ.
• Toppling — flexural toppling after Goodman & Bray: the joint dips into the face and (90° − joint dip) < (face dip − φ).

Joint orientations may be supplied from mapping, or estimated by clustering DTM surface normals over steep ground. That estimate assumes the exposed face is the joint plane — reasonable on a clean rock face, meaningless on talus or soil cover. On the reference site the recovered sets carry ±20–30° spread, which is what clustering a natural topographic surface yields rather than structural data.

Assumptions, collected

1. Rocks are point masses: no block shape, no rotation, no angular momentum, no fragmentation on impact.
2. Surface interaction is described by three global parameters (R_n, R_t, μ). There is no per-material map, so bedrock, talus and soil behave identically.
3. The terrain is rigid and unchanging during a run. The post-event surface is computed afterwards and never fed back.
4. The DTM is bare earth. Vegetation is neither an obstacle nor a damper.
5. Detachment is geometric: any cell steeper than a threshold can release. There is no stability criterion, no joint control on release, no triggering mechanism, and no temporal frequency.
6. Blocks are spheres for collision purposes, with equal masses in a swarm.
7. A barrier is a vertical plane of uniform height and uniform energy capacity, with no deformation, no post-impact geometry change, and no residual capacity after one breach.
8. Over a site-sized extent, a geographic CRS can be treated as a local tangent plane.
9. For change detection, the two surfaces differ by a rigid shift plus noise — no scale, rotation or non-rigid distortion is estimated.
10. For kinematics, joint sets are planar, persistent and continuous, with friction the only strength term (no cohesion, no water pressure).

Limitations, collected

1. Parameters are uncalibrated literature placeholders unless back-analysed against an observed event at the site.
2. Restitution is not identifiable from runout on terrain where rocks mostly roll — see above.
3. Results depend on grid resolution. Micro-topography below the cell size cannot influence a bounce, so runout on a coarse grid is smoother and generally slightly longer. Bounce height does not converge at all at these cell sizes.
4. Trajectory playback shows a sample (default 150 rocks), strided across the release area; the hazard rasters use every rock.
5. Rock markers are drawn inflated at wide zoom so they remain visible; overlap on screen is not overlap in the physics.
6. Collisions carry no spin and cannot fragment a block — they can only redirect and damp.
7. DoD volumes are a lower bound, biased low by an amount that depends on deposit geometry. Surveys must co-register to roughly 0.1 m for volumes to be trustworthy.
8. Hydrological flow terminates in filled flats (5% of cells on the reference DTM); D8 also over-concentrates flow on divergent hillslopes.
9. Kinematic "possible" is not "likely", and DTM-derived joint sets are a first look only.
10. One shared PIN, so results are attributable to the team, not to a person.
11. Structures are classified heuristically from the nDSM, and a bare-earth run (the default) over-predicts runout wherever the ground carries buildings or forest.
12. Flood screening is geometric, not hydraulic — no timing, momentum, backwater or structures, and rainfall mode's inputs are assumptions rather than measurements.
13. No slope stability, no factor of safety, no probability of failure is computed anywhere in this application.

Verification summary

Each module was checked against cases with a known answer rather than against plausibility.

Basis and references

The trajectory model follows the standard lumped-mass approach used by established rockfall codes (RocFall, Rockyfor3D). Default restitution and friction ranges are taken from Pfeiffer & Bowen (1989) and the Rockyfor3D / RocFall documented ranges: normal restitution ~0.40–0.53 on bedrock, ~0.30 on talus, ~0.25 on soil; tangential restitution ~0.85 on smooth rock, lower on rough or soft ground; rolling friction μ ≈ tan(29°).

Change detection follows the repeat-survey literature: Nuth & Kääb (2011) for co-registration, and the level-of-detection and volumetric error-budget approach of Brasington, Wheaton and co-workers. Sink filling uses the priority-flood algorithm of Barnes et al. (2014); D8 routing follows O'Callaghan & Mark (1984). Kinematic tests follow Markland, with the flexural toppling criterion of Goodman & Bray (1976).

Every default value is stated with its source in backend/params.py. Values are placeholders until site-specific data replaces them.

See the Guide for operation and the Changelog for what changed when — entries that alter results are marked.