DRAFT: Flexible threshold

[amoodie]

This was an earlier attempt at integrating some shear stress info to be able to remove the water and stop routing based on the stress.

Something is wrong though, because stress ends up being equal everywhere, so I obviously messed something up.

I think the first step would be to add another elseif branch here to do something like if grid.stess(indNew) < criticalstress; continuePropagateAvulsion = false; stop the routing. This would at least take care of the avulsion pathfinding stuff.

The next step would be to make the channel abandoning process work based on the stress.

@alimaye I didn't get very far in this buggy branch before trying to resolve other issues with avulsions, but I thought it might be helpful if you decide to take a crack at it while I'm in the field. You can either make a copy of my branch and then build on top of it and PR your own branch (the best-practices git workflow) or you could just start your own branch and then copy-paste the helpful bits.

[alimaye]

It may be simpler to differentiate the threshold discharge for abandoning a branch from the threshold stress for topographic change. Implementing the latter should be only require setting Qs = 0 for cells with grid.stess(indNew) < criticalstress

[amoodie]

So I'm still working on this, and have made some progress, but the critical stress thing won't work, I don't think. Revisitng the original paper, it's clear that the assumption is bankfull normal flow everywhere, where the width and depth are solved based on stress t*==1.8 everywhere.

I think the choice is then to use a critical flow depth to stop routing? But how to determine this depth...

[alimaye]

Ah that's right, the bankfull assumption is embedded. Could peg a critical flow depth to scale with the grain size D, which is specified as an input parameter. I don't know what to plug in for a scale factor but could make it a parameter for now.

[amoodie]

stashing the attempt I made to use a depth threshold, in case we want it for some reason in the future. I've decided to go another route, and will submit the PR shortly.

In propogate avulsions, one each step, compute the following estimate of depth

            dz = -(gridTemp.z(indNew)-gridTemp.z(indCurrent));
            [indCurrentRow,indCurrentCol] = ind2sub(grid.size,indCurrent); % current
            [indNewRow,indNewCol] = ind2sub(grid.size,indNew); % new
            dL = gridTemp.dx*sqrt((indCurrentRow-indNewRow)^2+(indCurrentCol-indNewCol)^2);
            gridTemp.S.alongFlow(indNew) = dz/dL;
            estimatedH = (alpha_b./gridTemp.S.alongFlow)*D; % depth

then check as a possible condition to stop routing

            % Stop path construction if the flow depth in the new cell
            % would be very small.
            elseif estimatedH(indNew) < 0.5

                % end iteration of the while loop
                continuePropagateAvulsion = false;