water_tools¶
The route_water
routine manages the water routing.
During route_water
, water iteration is repeated a total of itermax
times.
During each of these iterations of the water routing, the following methods are called in order:
Init the water iteration routine. 

Run a single iteration of travel paths for all water parcels. 

Calculate free surface after routing all water parcels. 


Finish updating flow fields. 
Public API methods attached to model¶
The following methods are defined in the water_tools
class, of the pyDeltaRCM.water_tools
module.
 class pyDeltaRCM.water_tools.water_tools¶
 build_weight_array(array, fix_edges=False, normalize=False)¶
Weighting array of neighbors.
Create np.array((8,L,W)) of quantity a in each of the neighbors to a cell
 check_for_boundary(inds)¶
Check whether parcels have reached the boundary.
Checks whether any parcels have reached the model boundaries. If they have, then update the information in
free_surf_flag
. Parameters
inds (
ndarray
) – Unraveled indicies of parcels.
 check_size_of_indices_matrix(it)¶
Check if step path matrix needs to be made larger.
Initial size of self.free_surf_walk_inds is half of self.stepmax because the number of iterations doesn’t go beyond that for many timesteps.
Once it reaches it > self.stepmax/2 once, make the size self.iter for all further timesteps
 compute_free_surface()¶
Calculate free surface after routing all water parcels.
This method uses the free_surf_walk_inds matrix accumulated during the routing of the water parcels (in
run_water_iteration
) to determine the free surface. The operations of the free surface computation are placed in a jitted function_accumulate_free_surface_walks
. Following this computation, the free surface is smoothed by steps infinalize_free_surface
.See 1 and 2 for a complete description of hydrodynamic assumptions in the DeltaRCM model.
Examples
The sequence of compute_free_surface is depicted in the figures below. The first image depicts the “input” to compute_free_surface, which is the current bed elevation, and the path of each water parcel (cyan lines in right image).
The compute_free_surface method then 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
).The output from
_accumulate_free_surface_walks
is then used to calculate a new stage surface (H_new
) based only on the water parcel paths and expected water surface elevations, approximately asHnew = sfc_sum / sfc_visit
(“computed Hnew” in figure below)Following this step, a correction is applied forcing the new free surface to be below sea level and above or equal to the land surface elevation over the model domain (“stage” in figure below). This surface Hnew is used in the following operation
finalize_free_surface
.Finally, the updated water surface is combined with the previous timestep’s water surface and an underrelaxation coefficient
omega_sfc
. 1
A reducedcomplexity 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/esurf3672015
 2
A reducedcomplexity 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/esurf3872015
 finalize_free_surface()¶
Finalize the water free surface.
This method occurs after the initial computation of the free surface, and creates the new free surface by smoothing the newly computed free surface with a jitted function (
_smooth_free_surface
), and then combining the old surface and new smoothed surface with underrelaxation. Finally, aflooding_correction
is applied.
 finalize_water_iteration(iteration)¶
Finish updating flow fields.
Clean up at end of water iteration
 flooding_correction()¶
Flood dry cells along the shore if necessary.
Check the neighbors of all dry cells. If any dry cells have wet neighbors, check that their stage is not higher than the bed elevation of the center cell. If it is, flood the dry cell.
 get_water_weight_array()¶
Get step direction weights for each cell.
This method is called once, before parcels are stepped, because the weights do not change during the stepping of parcels.
See Hydrodynamics for a description of the model design and equations/algorithms for how water weights are determined.
Note
No computation actually occurs in this method. Internally, the method calls the jitted
_get_water_weight_array
, which in turn loops through all of the cells in the model domain and determines the water weight for that cell withget_weight_sfc_int
and_get_weight_at_cell_water
.Examples
The following figure shows several examples of locations within the model domain, and the corresponding water routing weights determined for that location.
 get_wet_mask_nh()¶
Get wet mask.
Returns np.array((8,L,W)), for each neighbor around a cell with 1 if the neighbor is wet and 0 if dry
 init_water_iteration()¶
Init the water iteration routine.
 route_water()¶
Water routing main method.
This is the main method for water routing in the model. It is called once per update() call.
 Internally, this method calls:
 run_water_iteration()¶
Run a single iteration of travel paths for all water parcels.
Runs all water parcels (Np_water parcels) for stepmax steps, or until the parcels reach a boundary.
All parcels are processed in parallel, taking one step for each loop of the
while
loop.Example
On the initial delta surface, see how ten selected parcels are routed through the domain:
 update_Q(dist, current_inds, next_inds, astep, jstep, istep, update_current=False, update_next=False)¶
Update discharge field values.
Method is called after one set of water parcel steps.
 update_flow_field(iteration)¶
Update water discharge.
Update the water discharge field after one set of water parcels iteration.
 update_velocity_field()¶
Update flow velocity fields.
Update the flow velocity fields after one set of water parcels iteration.
water_tools helper functions¶
The following routines are jitted for speed. They generally take a large number of arrays and constants and return a new array(s) to continue with the model progression within the methods defined above.
 pyDeltaRCM.water_tools._get_weight_at_cell_water(ind, weight_sfc, weight_int, depth_nbrs, ct_nbrs, dry_depth, gamma, theta)¶
Compute water weights for a given cell.
This is a jitted function called by
_get_water_weight_array()
.
 pyDeltaRCM.water_tools._choose_next_directions(inds, water_weights)¶
Get new cell locations, based on water weights.
 Algorithm is to:
loop through each parcel, which is described by a pair in the inds array.
determine the water weights for that location (from water_weights)
choose a new cell based on the probabilities of the weights (using the random_pick function)
 Parameters
inds (
ndarray
) – Current unraveled indices of the parcels.(N,)
ndarray containing the unraveled indices.water_weights (
ndarray
) – Weights of every water cell.(LxW, 9)
ndarray, uses unraveled indicies along 0th dimension; 9 cells represent self and 8 neighboring cells.
 Returns
next_direction – The direction to move towards the new cell for water parcels, relative to the current location. I.e., this is the D8 direction the parcel is going to travel in the next stage,
pyDeltaRCM.shared_tools._calculate_new_ind
. Return type
ndarray
 pyDeltaRCM.water_tools._calculate_new_inds(current_inds, new_direction, ravel_walk)¶
Calculate the new location (current_inds) of parcels.
Use the information of the current parcel (current_inds) in conjunction with the D8 direction the parcel needs to travel (new_direction) to determine the new current_inds of each parcel.
In implementation, we use the flattened ravel_walk array, but the result is identical to unraveling the index, adding iwalk and jwalk to the location and then raveling the index back.
ind_tuple = shared_tools.custom_unravel(ind, domain_shape) new_ind = (ind_tuple[0] + jwalk[newd], ind_tuple[1] + iwalk[newd]) new_inds[p] = shared_tools.custom_ravel(new_ind, domain_shape)
 pyDeltaRCM.water_tools._check_for_loops(free_surf_walk_inds, new_inds, _step, L0, CTR, stage_above_SL)¶
Check for loops in water parcel pathways.
Look for looping random walks. I.e., this function checks for where a parcel will return on its
new_inds
to somewhere it has already been infree_surf_walk_inds
. If the loop is found, the parcel is relocated along the mean transport vector of the parcel, which is computed as the vector from the cell (0, CTR) to the new location in new_inds.This implementation of loop checking will relocate any parcel that has looped, but only disqualifies a parcel p from contributing to the free surface in
_accumulate_free_surf_walks
(i.e., looped[p] == 1) if the stage at the looped location is above the sea level in the domain. Parameters
free_surf_walk_inds – Array recording the walk of parcels. Shape is (:obj:`Np_water, …)`, where the second dimension will depend on the step number, but records each step of the parcel. Each element in the array records the flat index into the domain.
new_inds – Array recording the new index for each water parcel, if the step is taken. Shape is (Np_water, 1), with each element recording the flat index into the domain shape.
_step – Step number of water parcels.
L0 – Domain shape parameter, number of cells inlet length.
CTR – Domain shape parameter, index along inlet wall making the center of the domain. I.e., (0, CTR) is the midpoint across the inlet, along the inlet domain edge.
stage_above_SL – Water surface elevation minuns the domain sea level.
 Returns
new_inds – An updated array of parcel indicies, where the index of a parcel has been changed, if and only if, that parcel was looped.
looped – A binary integer array indicating whether a parcel was determined to have been looped, and should be disqualified from the free surface computation.
Examples
The following shows an example of how water parcels that looped along their paths would be relocated. Note than in this example, the parcels are artificially forced to loop, just for the sake of demonstration.
 pyDeltaRCM.water_tools._update_dirQfield(qfield, dist, inds, astep, dirstep)¶
Update unit vector of water flux in x or y.
 pyDeltaRCM.water_tools._update_absQfield(qfield, dist, inds, astep, Qp_water, dx)¶
Update norm of water flux vector.
 pyDeltaRCM.water_tools._accumulate_free_surface_walks(free_surf_walk_inds, free_surf_flag, cell_type, uw, ux, uy, depth, dx, u0, h0, H_SL, S0)¶
Accumulate the free surface by walking parcel paths.
This routine comprises the hydrodynamic physicsbased computations.
 Algorithm is to:
loop through every parcel’s directed random walk in series.
for a parcel’s walk, unravel the indices and determine whether the parcel should contribute to the free surface. Parcels are considered contributors if they have reached the ocean and if they are not looped pathways.
then, we begin at the downstream end of the parcel’s walk and iterate upwalk until, determining the Hnew for each location. Downstream of the shorelineocean boundary, the water surface elevation is set to the sea level. Upstream of the shorelineocean boundary, the water surface is determined according to the landslope (
S0
) and the parcel pathway.repeat from 2, for each parcel.
Examples
The following shows an example of the walk of a few water parcels, along with the resultant computed water surface.
 pyDeltaRCM.water_tools._smooth_free_surface(Hin, cell_type, Nsmooth, Csmooth)¶
Smooth the free surface.
 Parameters
Hin – Stage input to the smoothing (i.e., the old stage).
cell_type – Type of each cell.
Nsmooth – Number of times to smooth the free surface (i.e., stage)
Csmooth – Underrelaxation coefficient for smoothing iterations.
Examples
The following shows an example of the water surface smoothing process.