pyDeltaRCM approximates hydrodynamics through the use of a weighted random walk. See 1 and 2 for a complete description of hydrodynamic assumptions in the DeltaRCM model. In this documentation, we focus on the details of model implementation, rather than model design.
Routing individual water parcels¶
Probabilities for water parcel routing to all neighbors and to self for each cell are computed once at the beginning of the water routing routine (
get_water_weight_array called from
Water routing probability for a given cell to neighbor cell is computed according to:
where is the local routing direction and is a unit vector pointing to neighbor from cell , and is the cellular distance to neighbor ( for cells in main compass directions and for corner cells. is a flow resistance estimated as an inverse function of local water depth ():
The exponent takes a value of unity by default for water routing (
theta_water, Default Model Variable Values), leading to routing weight for neighbor cell :
These weights above are calculated only for wet neighbor cells; all dry neighbor cells take a weight value of 0 (
Finally, probability for routing from cell to cell is calculated as:
Weights are accumulated for 8 neighbors and a probability of 0 is assigned to moving from cell to cell (i.e., no movement).
These 9 probabilities are organized into an array
self.water_weights with shape (
The following figure shows several examples of locations within the model domain, and the corresponding water routing weights determined for that location.
Because probabilities are computed for all locations once at the beginning of water iteration, all water parcels can be routed in parallel step-by-step in
During iteration, the direction of the random walk is chosen for each parcel via
_choose_next_directions, which internally uses
random_pick() for randomness.
For example, see the random walks of several parcels below:
add sentence or two above about check_for_loops.
Water routing completes when all water parcels have either 1) reached the model domain boundary, 2) taken a number of steps exceeding
stepmax, or 3) been removed from further routing via the
Combining parcels into free surface¶
Following the routing of water parcels, these walks must be converted in some meaningful way to a model field representing a free surface (i.e., the water stage).
compute_free_surface() is called, which takes as input the current bed elevation, and the path of each water parcel (top row in figure below).
compute_free_surface() method internally calls the
_accumulate_free_surface_walks() function to determine 1) the number of times each cell has been visited by a water parcel
sfc_visit), and 2) the total sum of expected elevations of the water surface at each cell (
_accumulate_free_surface_walks() itself iterates through each water parcel, beginning from the end-point of the path, and working upstream; note that parcels that have been determined to “loop” (
_check_for_loops() and described above) are excluded from computation in determining the free surface.
While downstream of the land-ocean boundary (determined by a depth-or-velocity threshold), the water surface elevation is assumed to be
0, whereas upstream of this boundary, the predicted elevation of the water surface is determined by the distance from the previously identified water surface elevation and the background land slope (
S0), such the the water surface maintains an approximately constant slope for each parcel pathway.
The algorithm tracks the number of times each cell has been visited by a water parcel (
sfc_visit), and the total sum of expected elevations of the water surface at each cell (
sfc_sum), by adding the predicted surface elevation of each parcel step while iterating through each step of each parcel.
Next, the output from
_accumulate_free_surface_walks() is used to calculate a new stage surface (
H_new) based only on the water parcel paths and expected water surface elevations, approximately as
H_new = sfc_sum / sfc_visit.
The updated water surface is combined with the previous timestep’s water surface and an under-relaxation coefficient (
With a new free surface computed, a few final operations prepare the surface for boundary condition updates and eventually being passed to the sediment routing operations.
A non-linear smoothing operation is applied to the free surface, whereby wet cells are iteratively averaged with neighboring wet cells to yield an overall smoother surface.
The smoothing is handled by
_smooth_free_surface() and depends on the number of iterations (
Nsmooth) and a weighting coefficient (
flooding_correction() is applied to the domain.
In this correction, all “dry” cells (a cell where the flow depth is less than the dry_depth) are checked for any neighboring cells where the water surface elevation (stage) is higher than the bed elevation of the dry cell.
If this condition is met for a given dry cell, the dry cell is flooded: the stage of the dry cell is set to the maximum stage of neighboring cells.
Finalizing and boundary conditions to sediment routing¶
Incomplete. Need to describe the updating of depth from stage, limiting everything to above H_SL, updating velocity and discharge fields, etc.
A reduced-complexity model for river delta formation – Part 1: Modeling deltas with channel dynamics, M. Liang, V. R. Voller, and C. Paola, Earth Surf. Dynam., 3, 67–86, 2015. https://doi.org/10.5194/esurf-3-67-2015
A reduced-complexity model for river delta formation – Part 2: Assessment of the flow routing scheme, M. Liang, N. Geleynse, D. A. Edmonds, and P. Passalacqua, Earth Surf. Dynam., 3, 87–104, 2015. https://doi.org/10.5194/esurf-3-87-2015