Estimating the volume of the South Lhonak Lake outburst

The Oct 4, 2023, flashflood in Sikkim was caused by the outburst of South Lhonak Lake. In my previous post, I showed how to derive some basic metrics about the outburst path and its length to the destroyed Chungthang dam which is an integral part of the Teesta III 1200 Mw hydropower project. In the meantime, Planet satellite imagery before and after the event offer insights into what happened at the lake itself. A number of evidences propose that the lake may have burst due to a failure of the moraine straddling the lake in the northern part.

To simulate such floods requires good estimates of the outflow volume and while we are still waiting for measurements taken in the field, I will try to come up with an estimate based on topographic analysis of the High Mountain Asia (HMA) 8-m DEM together with Planet imagery taken before and after the disaster. The following post will detail the steps that I have taken in GIS and TopoToolbox, and will explain the assumptions behind the analysis. For those who do not want to read the whole thing, here’s the result of the analysis: The outburst volume is 41.4 Mio m³, which equals 16544 Olympic size swimming pools.

Ok, so I started with three different datasets which I used to estimate the volume and which are shown below:

1. The HMA DEM
2. Planet imagery of the lake before the outburst
3. Planet imagery of the lake after the outburst

In QGIS, I performed following steps.

1. Loading the DEM.
2. Reprojecting the DEM so that it has the same projection like the Planet imagery
3. Digitizing the outlines of the lakes before and after the outburst (of course, automated detection of water surfaces could have also been possible, but a lot of ice on the drained lake may have complicated that step)
4. Export all data so that I can work with it in TopoToolbox.

In TopoToolbox, I loaded and preprocessed the data using following code (I am not going to explain each step in detail, but comment the code below).

DEM = GRIDobj('HMA_DEM_utm45.tif');
% There are some missing data which we'll fill with interpolation
DEM = inpaintnans(DEM);

lake       = shaperead("lake.shp");
lake_pre   = polygon2GRIDobj(DEM,lake(1));
lake_post  = polygon2GRIDobj(DEM,lake(2));

The DEM is derived from optical imagery that is a couple of years older than the pre-image. To remove the ice-elevations that are now ice-free and waterfilled, I assigned the lowest DEM value to all DEM pixels contained in the lake_pre mask.

DEM.Z(lake_pre.Z) = min(DEM.Z(lake_pre.Z));

The image below shows the original DEM and the DEM adjusted for the pre-outburst situation.

Now, how can one estimate the elevation of the lake after the outburst? The moraine along the southern part of the lake has a very straight slope and here I simply assumed that the slope can be linearly extrapolated to the shoreline following lake drainage. I drew a couple of start and end points to draw straight lines along which I estimate the slope. Here are the linear indices of the start and end points which I then convert to pairs of coordinates.

ix = [389062
      385693
      453271
      443878
      469983
      458582
      484688
      471279
      429865
      423148
      405785
      401071];

[x,y] = ind2coord(DEM,ix);

x = reshape(x,2,[])';
y = reshape(y,2,[])';

The next figure shows the DEM with the profiles.

Then, I created a new DEM with the lake surface removed:

DEMc = clip(DEM,~lake_pre);

Now, here is the part of the code, in which I take the start and end points of the profiles, extract elevation values from the clipped DEM, and then fit a linear model to each line:

figure
ax = axes;
for r = 1:size(x,1)
    p = demprofile(DEMc,50,x(r,:),y(r,:));
    I = ~isnan(p.z);
    mdl = fitlm(p.d(I),p.z(I));

    % find distance of the contour line of the drained lake
    [~,~,II] = polyxpoly(p.x,p.y,lakes(2).X,lakes(2).Y); 
    
    di = p.d(II(1));
    [zpred(r),zpredci(r,:)] = predict(mdl,di);
    
    hold on
    plot(p.d,p.z,'-k')
    errorbar(di,zpred(r),zpred(r)-zpredci(r,1),zpredci(r,2)-zpred(r));
     
end

box on
yline(mean(zpred),'b--','Lakelevel (post outburst)',...
    'LabelHorizontalAlignment','left')
xlabel('Distance along profiles [m]')
ylabel('Elevation [m]')

zpred and zpredci (I won’t use those here) contain the estimated elevations of the shoreline of the lake after drainage.

Unfortunately, the lake is difficult to map close to the glacier terminus, hence I am using a little trick to derive the elevations in these areas (based on the distance transform inward into the lake).

D = distance(DEM,~lake_pre);
[X,~] = getcoordinates(D);
DEMc.Z(lake_post.Z | (D.Z>100 & X < 617500)) = mean(zpred);
DEMc.Z(bwperim(lake_pre.Z)) = DEM.Z(bwperim(lake_pre.Z)); 

Now we have an estimate of the elevation for the lake area after drainage, and we can use Laplacian Interpolation to estimate elevations between the before- and after-shorelines.

DEMc = inpaintnans(DEMc);

The post-DEM then looks as below (with pre-DEM for comparison):

Based on these DEMs, we can calculate a DEM of difference and calculate the volume lost.

nansum(DEM.Z(:)-DEMc.Z(:))*DEM.cellsize^2 
ans =

   41361412

The volumetric difference is 41.4 Mio m³, which equals to 16544 Olympic size swimming pools.

This volume, however, does not include the eroded volume, in particular along the spill way. We can artificially cut the DEM using minima imposition.

FD = FLOWobj(DEMc);
DEMcc = imposemin(FD,DEMc);
imageschs(DEMcc,[],'colormap','landcolor')

To wrap up, this analysis relied on the assumption that we can linearly extrapolate the straight slopes on the southern bank of the lake to estimate the lake level after the outburst. The part of the lake close to the terminus is difficult to reconstruct since ice in lake after the outburst complicates deriving a good outline. To this end, the calculations yield an estimate of the outburst volume of 41.4 Mio cubicmeters. This point estimate can be enhanced by confidence intervals when we account for the variability of the mean of estimates of lake-level elevations. But this might be a topic for an upcoming blog post.