Implement avalanche2d

[wkearn]

For landscape evolution modeling purposes, we would like to implement the following MATLAB function avalanche2d in libtopotoolbox. The idea is to redistribute loose sediment over a fixed bedrock surface so that the angle of repose constraint is satisfied. This should be fairly straightforward to port in its current form, though I do wonder whether there are any algorithmic improvements to be made, in particular, whether one could limit the number of iterations in a way similar to how we use the fast marching method in excesstopography_fmm2d.

We should probably start with a signature like

TOPOTOOLBOX_API
void avalanche2D(float *out_thickness, float *dem, float *in_thickness, 
                            float threshold_slope, float minimum_thickness,
                            float cellsize, ptrdiff_t nrows, ptrdiff_t ncols);

We would also like to handle boundary conditions somewhat more flexibly, to assign Neumann or Dirichlet BCs on different segments of the DEM boundary. This will require some thought about the interface, so it may be something to consider after we have the initial implementation working.

function [H] = avalanche2d(Z,H,dx,Sc,minH)
% function to move loose material (e.g., snow) downhill and deposit at
% angle of repose
% Adapted from Mark Kessler's GC2D program
% <Z> is the elevation of stable material, e.g., the bedrock surface
% <H> is the thickness of the mobile material, e.g., the snow depth
% <Sc> is the angle of repose

% Calculate elevation of loose material
Zi = Z + H;

% Adjust angle of repose (from degrees to m/m)
dHRepose = dx * tan(deg2rad(Sc)) ;

% Size of model domain
[rws,cls] = size(Z) ;

Ho = max( H, 0 ) ; % Negative thicknesses not allowed

% Count number of iterations
ct = 1;

while ct<1e2 %(1)
    % calculate elevation drop of loose material surface in x direction
    dZidx_down = zeros(rws,cls);
    dZidx_down(:,2:end) = max( 0, ( Zi(:,2:end) - Zi(:,1:end-1) ) ) ;
    dZidx_up = zeros(rws,cls);
    dZidx_up(:,1:end-1) = max( 0, ( Zi(:,1:end-1) - Zi(:,2:end) ) ) ;
    dZidx = max( dZidx_up, dZidx_down ) ;
    % calculate elevation drop of loose material surface in y direction
    dZidy_left = zeros(rws,cls);
    dZidy_left(2:end,:) = max( 0, ( Zi(2:end,:) - Zi(1:end-1,:) ) ) ;
    dZidy_right = zeros(rws,cls);
    dZidy_right(1:end-1,:) = max( 0, ( Zi(1:end-1,:) - Zi(2:end,:) ) ) ;
    dZidy = max( dZidy_left, dZidy_right ) ;
    % calculate mean elevation drop of loose material surface
    grad = sqrt( dZidx.^2 + dZidy.^2 );
    gradT =  dZidy_left + dZidy_right + dZidx_down + dZidx_up ;
    gradT(gradT==0) = 1;
    grad(Ho<minH) = 0 ;

    % Exit condition
    mxGrad = max( max( grad ) ) ;
    if ( mxGrad <= 1.1*dHRepose )
        break ; % avalanching condition not given (anymore)
    end

    % ???
    delH = max( 0, ( grad - dHRepose ) / 3 ) ; 

    % ???
    Htmp = Ho ;
    Ho = max( 0, Htmp - delH );
    delH = Htmp - Ho ;

    % ???
    delHdn = zeros(rws,cls) ; delHup = zeros(rws,cls) ;
    delHlt = zeros(rws,cls) ; delHrt = zeros(rws,cls) ;

    delHup(:,2:end) = delH(:,1:end-1) .* dZidx_up(:,1:end-1)./gradT(:,1:end-1) ;
    delHdn(:,1:end-1) = delH(:,2:end) .* dZidx_down(:,2:end)./gradT(:,2:end) ;
    delHrt(2:end,:) = delH(1:end-1,:) .* dZidy_right(1:end-1,:)./gradT(1:end-1,:) ;
    delHlt(1:end-1,:) = delH(2:end,:) .* dZidy_left(2:end,:)./gradT(2:end,:) ;

    Ho = max( 0, Ho + delHdn + delHup + delHlt + delHrt ) ;

    Zi = Z + Ho ;

    ct = ct+1;
end

H = Ho + (H<0).*H ;

end % of function