Hydrodynamics¶
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.
The water routing operations are orchestrated by route_water; the following sections narrate the sequence of events within a single call of this method.
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 run_water_iteration).
Water routing probability for a given cell \(j\) to neighbor cell \(i\) is computed according to:
where \(\mathbf{F}\) is the local routing direction and \(\mathbf{d_i}\) is a unit vector pointing to neighbor \(i\) from cell \(j\), and \(\Delta_i\) is the cellular distance to neighbor \(i\) (\(1\) for cells in main compass directions and \(\sqrt{2}\) for corner cells. \(R_i\) is a flow resistance estimated as an inverse function of local water depth (\(h_i\)):
The exponent \(\theta\) takes a value of unity by default for water routing (theta_water, Default Model Variable Values), leading to routing weight for neighbor cell \(i\):
These weights above are calculated only for wet neighbor cells; all dry neighbor cells take a weight value of 0 (_get_weight_at_cell_water).
Finally, probability for routing from cell \(j\) to cell \(i\) is calculated as:
Weights are accumulated for 8 neighbors and a probability of 0 is assigned to moving from cell \(j\) to cell \(j\) (i.e., no movement).
These 9 probabilities are organized into an array self.water_weights with shape (L, W, 9)`.
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 run_water_iteration.
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:
Todo
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 _check_for_loops function.
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).
First, the 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).
The 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 (sfc_sum).
_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 (omega_sfc).
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 (inside finalize_free_surface()).
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 (Csmooth).
Finally, a 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.
Similar to _smooth_free_surface() described above, flooding_correction() acts to remove roughness in the water surface and contributes to model stability.
Finalizing and boundary conditions to sediment routing¶
The final step of the water routing operations in each call to route_water is the updating of model fields, and application of a few corrections and boundary conditions.
These operations are handled within finalize_water_iteration.
First, the stage field is limited to elevations above or equal to sea level H_SL. Then, the depth field is updated as the vertical distance between the stage field and the bed elevation eta. These operations ensure that these fields are in agreement over the entire model domain.
Next, methods update_flow_field and update_velocity_field are called, which handle the updating of water discharge and water velocity fields, respectively.
Note
The next step in the model update sequence is sediment routing (route_sediment). For more information on the next stage, see Morphodynamics.