Remote Sensing-Derived Bathymetry of Lake Poopó

Located within the Altiplano at 3,686 m above sea level, Lake Poopó is remarkably shallow and very sensitive to hydrologic recharge. Progressive drying has been observed in the entire Titicaca-Poopó-Desaguadero-Salar de Coipasa (TPDS) system during the last decade, causing dramatic changes to Lake Poopó’s surface and its regional water supplies. Our research aims to improve understanding of Lake Poopó water storage capacity. Thus, we propose a new method based on freely available remote sensing data to reproduce Lake Poopó bathymetry. Laser ranging altimeter ICESat (Ice, Cloud, and land Elevation Satellite) is used during the lake’s lowest stages to measure vertical heights with high precision over dry land. These heights are used to estimate elevations of water contours obtained with Landsat imagery. Contour points with assigned elevation are filtered and grouped in a points cloud. Mesh gridding and interpolation function are then applied to construct 3D bathymetry. Complementary analysis of Moderate Resolution Imaging Spectroradiometer (MODIS) surfaces from 2000 to 2012 combined with bathymetry gives water levels and storage evolution every 8 days.


Introduction
In the present study, we propose a new method of bathymetry restitution developed for the Lake Poopó in Bolivia.With this method we have recovered Lake Poopó water levels and storage time series which have remained unknown until now.The method is easy to implement and is based on freely available remote sensing data.Geographical applicability of our method is however limited to lakes situated in arid or semi-arid regions.Lake Poopó's system (17°10'-19°53'S and 67°50'-66°24'W, at 3,686 m above sea level) is comprised of two interconnected lakes, Lake Uru-Uru and Lake Poopó (Figure 1a).Lake Poopó is remarkably shallow and sensitive to climate variations.Its shallow zone is seasonally dependent on hydrologic recharge during the wet season, which varies in extent from almost 0 to 3,000 km 2 [1].Lake Poopó is part of the Titicaca-Poopó-Desaguadero-Salar de Coipasa water system (TPDS) and plays an important role in the regional water budget.Unfortunately, Lake Poopó water levels and storage are not recorded.Previous works have attempted to correlate the dimensions of Lake Poopó with its depth in a relatively linear fashion [2,3], or have estimated Lake Poopó levels directly from Lake Titicaca level variations [2,4].Accurate information about the water stage or water supplies was not available in these studies.Benefiting from years of experience, space borne altimetry offers the best solution for ungauged or difficult to access lakes and rivers.The usual approach consists in establishing remote sensed hypsometry of a lake, based on the simple assumption that every water level change induces a change of the water surface.If the water surface variation is large enough, it can be easily detected with multispectral high or medium resolution imagery [5][6][7][8].When correlated with level variation time series or topographic measurements, the surface information can be converted into volume [3,6,9].This technique was successfully applied in some studies [5,6] but unfortunately, Topex/Poseidon (T/P) is the only reliable radar altimeter satellite to cover Lake Poopó, and then only on its shallow western side, which is intermittently covered by water.While topographic, sonar and airborne lidar measurements are costly and difficult to obtain [3] it remains impossible to produce regular time series of Lake Poopó levels with radar altimetry.Nevertheless, several other methods and instruments can be used to create a lake bathymetry.One of these methods is very high resolution stereo imaging [10] but this method cannot be applied over flat and uniform surfaces of the lake bed.Photobathymetry [11,12] can be a solution for very shallow lakes.This method requires high resolution imaging and several assumptions about water turbidity and bottom reflectance.In the case of Lake Poopó these limitations can lead to a bathymetry with non-uniform and average vertical precision.We also examined the possibility of applying SAR imaging.However, this technique does not provide enough information in terms of vertical precision [13].In looking for centimetric precision, we propose using the first laser-ranging instrument GLAS onboard the ICESat mission.ICESat provides even better accuracy than radar altimetry in non-clouded conditions [14][15][16][17] and has already been applied in lake surface/depth studies [18,19].A rapid examination of ICESat tracks shows their perfect location just in the middle of Lake Poopó (Figure 1b).With a 91-day orbit it is possible to produce a limited time series between 2004 and 2009.In [16], the ICESat data have been validated over the Salar de Uyuni and it was found that the accuracy of laser ranging measurements is ±3 cm over dry land.Our research aims to exploit this laser-ranging high precision to estimate elevations of water contours obtained with Landsat imagery and transform them into isolines (isobaths).In another study, Abileah and Vignudelli [20] propose the combination of radar altimetry and imagery to calculate Lake Nasser levels and depth contours.In contrast, our method takes laser altimetry measurements directly over the bare ground when the lake reaches a very low stage.This approach avoids uncertainties due to wave height or wind direction because laser altimetry measures vertical height with high precision over dry land.Finally complementary analysis of Moderate Resolution Imaging Spectroradiometer (MODIS) surfaces from 2000 to 2012 combined with bathymetry gives water levels and storage evolution every 8-day.The article is organized as follows: Section 2 describes the data Section 3 describes the methodology for estimating contour elevations.Section 4 presents the results.Conclusions are drawn in Section 5.

Data
Through the USGS GLOVIS images archive [21] 18 cloud free scenes from Landsat 7 ETM+ and Landsat 5 TM were selected between 1999 and 2008 (Section 3.2).All selected Landsat scenes have Standard Terrain Correction (Level 1T).Image projection is UTM zone 19N, with 30 × 30 m of ground resolution.Since the first Landsat mission in 1972, there are over 150 suitable scenes now available at no cost through the Landsat archives.Terra Moderate Resolution Imaging Spectroradiometer (MODIS) was launched in 1999 on the Terra satellite.Between 2000 and 2012, 560 MODIS scenes were downloaded [22].Data from the first laser-ranging instrument GLAS on board the ICESat 2003-2009 mission is used [23].With a 91-day orbit, the ICESat GLA06 v33 [17] global elevation product is used, prior to the Correction to ICESat Data Product Surface Elevations announced in January 2013 [24].It provides laser footprint geolocation and surface elevations determined with ice sheet specific algorithms [25].Elevation values are referred to the T/P ellipsoid.GLA14 dataset showing only 1/3 of available flyovers was discarded.

General Introduction
In this study, we use laser ranging information to convert water masks obtained with Landsat imagery (Section 3.2) into isolines (Section 3.4).Isolines (isobaths) are obtained by extracting water area contours and by interpolating ground elevation values in each location where water contours intersect with ICESat elevation profiles.Individual points from isolines are grouped in a cloud of 3D points and a bathymetry is created using grid meshing and interpolating function (Section 3.5).The bathymetry was created in the limits of existing ICESat profiles and depths below 3,685.70m are not shown.Moreover, between 2000 and 2012, complementary analysis of variations of water area inferred from MODIS images classification (Section 3.2) was performed and linked to bathymetry through the Landsat/MODIS correlation curve.Lake Poopó water levels and storage are retrieved by combining bathymetry and MODIS surfaces time series.

Water Surface Masks for Bathimetry and Water Surface Area Variations
The Landsat 5 and Landsat 7 images selected for this study were strictly limited to cloud free and inundation free scenes in order to reduce difficulties in processing and interpretation.In order to extract water surfaces from Landsat images, digital numbers in each spectral band were converted to surface reflectance using the calibration coefficients provided in the image's metadata.The water masks were obtained by combining supervised classification of the Normalized Difference Water Index (NDWI) [26], Modified Normalized Difference Water Index (MNDWI) [27], Weighted Water Index [28,29] and single band thresholding [30,31].Since the threshold value of the NDWI for delineating open water features is known to vary in multi-temporal studies [27,32] we manually adjusted the NDWI threshold on a satellite-by-satellite basis.Manual adjustment of the threshold could achieve a more accurate result in the water delineation [27,33,34].Co-registered water masks were sorted by decreasing surface order and then superposed.Any incoherencies between consecutive masks, e.g., when a lower surface contour exceeded a greater surface contour, were binary eroded to leave at least one pixel of spacing between masks.Contours from different images were converted into two datasets: one representing the main water body contours, the other representing major interior details such as islands or sandbars.Moreover, the MODIS instrument provides new remote sensing opportunities for global daily monitoring with improved spatial resolution for land observation at 250 and 500 m and improved spectral band placement [35].The MODIS products include seven bands in the visible to middle infrared.Shallow depths such as those observed over the Lake Poopó; increase considerably the amount of solar energy reflected by a water body [36,37] showed that the strong water absorption at wavelength > 1 μm in MODIS bands (number 5, 6 and 7) does not allow illumination of the sediments in the water or at the shallow bottom of a water column.In this work a simple unsupervised classification based on threshold technique was performed on MODIS Band 5 to delineate the shallow, open water of the Lake Poopó.It has been assumed that for small value of surface reflectance allows characterizing open water, without any consideration on any normalized water index.This method was developed for floodplains in Arid zones [7] and has been tested and applied in other studies [1,[38][39][40][41][42]. Surface analysis of a set of 550 8-day MODIS mosaics was strictly limited to the main water body.Flood plains and rivers outlets were ignored.Finally, as MODIS and Landsat images have different resolutions, we established link between MODIS surfaces and our bathymetry through Landsat/MODIS correlation curve in Section 4. Surfaces lower than 450 km 2 are not used in this study, and are not shown on correlation and hypsometry curves.

ICESat Profiles Selection
In Figure 2, typical ICESat profiles from the northern part of Lake Poopó are shown.Water surface can be identified in the lowest part of each profile.Several outliners' points with differences of 20 to 40 cm to the water level can be attributed to waveform saturation [43].Further examination coupled with optical imagery permit to localize all saturated points in water-ground/mud transition zones (Figure 3).Saturation elevation correction (i_satElevCorr) as suggested by the GLAS Altimetry Data Dictionary [44] was applied.However, in transition zones this correction can produce sporadically underestimated elevations.All elevation profiles taken during or immediately after the February-March rain/flood season were rejected.These profiles can be easily identified by the presence of signal saturation zones (Figure 2b).The water stage of each ICESat profile is represented by a threshold value that lays 5 cm above the water surface (Figure 2a,b and Figure 3).This approach allows handling saturation problems occurring inside the main water body at the intersection of islands, dunes and capes.During processing, this threshold value is used to exclude from interpolation (in Section 3.4) all points that lay at the water level (points below threshold).Last, for all retained ICESat profiles (Table 1), elevation values referred to the T/P ellipsoid are converted to EGM2008 geoidal heights.

Creating Isolines (Isobaths)
For each main water body contour, a simple algebraic method of line segment intersection [45,46] is used to determine geographic coordinates where water contour and ICESat elevation profiles intersect.For each intersection several elevation values are linearly interpolated using all ICESat profiles.If (for a given ICESat profile), at least one measurement lies above the "water" threshold, the interpolated value is kept.Using this approach, at least several elevation values can be obtained for each main water body contour.All elevation values obtained within one intersection were filtered to discard erroneous points.Lake Poopó is exceptionally shallow and given large values of possible outliers in relation to the small number of valid points, the median offers a more robust predictor than the mean.Thus the median value is retained as final isoline elevation, and is also assigned to all interior contours of the same layer.All isolines, except the bottom one, were obtained with previously described methodology.For the lowest stage of the lake we first identified 5 October 2009 ICESat flyover as the lowest valid elevation profile.We then extracted a water surface mask from Landsat 7 image taken on 2 October 2009 and we completed it with Landsat 5 image from 10 October 2009.The elevation for the lowest stage mask was obtained by calculating the median value of all ICESat elevation points lying within the body of water.

Bathymetry
The bathymetry reconstruction was performed with Matlab functions: TriScatteredInterp function [47], allowing interpolation on a regular grid from irregularly scattered points, and MeshGrid function for grid meshing [48].A 5 × 5 pixel median filter was applied to smooth irregularities and remove artifacts.The bathymetry keeps the same pixel size as input Landsat images.For the purposes of comparison, Lake Poopó hypsometry is estimated directly from bathymetry using Heron's formula [8].The methodology presented above allowed us to calculate the Lake Poopo bathymetry and water storage variations as shown in the following section.

Results
Lake Poopó hypsometry is presented on Figure 4.The lowest surface attributed to the 3,685.70m elevation (ICESat 5 October 2009 flyover) is 467 km 2 while the highest point on the curve has 3,688.10m of elevation for a surface of 2,316 km 2 .The MODIS surfaces (Figure 5) are linked to bathymetry through the Landsat/MODIS calibration curve (Figure 6) and a correlation coefficient of R 2 = 0.97 was obtained.A top and inclined view of the final bathymetry is presented in Figure 7a,b.Table 2 shows relations between isolines vertical elevation, the elevation standard deviation (SD) and the number of profiles used for the calculation.The SD increases with the number of profiles and reaches 0.11 m for the highest contour.Vertical SD of the bottom mask calculated separately is 0.05 m.The precision of our method and repeatability has been confirmed by processing masks with similar water surface extension from images taken at different dates (a1-a2, b1-b2 in Table 2).In both cases, identical elevations were obtained.The bathymetry vertical resolution can be considered as decimetric.The bathymetry combined with MODIS time series is used to calculate Lake Poopó water levels and water storage (Figure 8).We can observe that lake maximum water storage is always reached between 1 March and 20 March each year.In March 2001 we observed the highest relative water storage.Minimum water storage can be observed at the beginning of November if only the water level in March is higher than 3,686.6 m (surface larger than 1,490 km 2 ).When the water level in March was lower than 3,686.6 m then the minimum lake storage was reached during the 1st week of October as it has been observed in 2009, 2010 and 2011.
From MODIS imagery (Figure 5) we can deduce that in 2009, 2010 and 2011 the water storage in October/November was close to zero.The bathymetry was created in the limits of existing ICESat profiles and depths below 3,685.70m are not shown.We estimate that the missing part of the storage represents 10%-15% of the total volume.

Discussion and Perspectives
Our method uses exclusively freely available remote sensing data, and the vertical standard deviation of our estimates is ±10 cm.All major interior details of the lake, such as islands and sandbars, can be reproduced if needed.The method presented here is somewhat limited.Availability of appropriate images with visible water surface changes and accurate contour detection are prerequisite.The presence of an elevation profile that coincides with the lake minimum stage is the most important requirement.In Lake Poopó's case the lowest ICESat profile was taken in October 2009.Later in November-December same year the lake dried out entirely (Figure 5).As for the method itself, the interpolation of elevation values in some regions of a lake can give incoherent results (sloped but flooded area seen from nadir, river path) and should be carefully analyzed.In our first attempt, we tried to split points into different zones (Figure 1b) and recommend this approach for future studies.In each zone, it is possible to have multiple intersections.For each contour, median values from different zones can be compared.Intersections on a sloping flooded area seen from nadir or in the river path can give values 1-3 m above the real water elevation.Those points should be rejected.Elevation values in opposite zones, i.e., south-north and east-west need to be checked.If the elevation values in opposite zones differ systematically from layer to layer in the same direction, this can be interpreted as a local geoid effect.If elevation values in opposite zones differ occasionally, the main water body contour may be affected by wind.In this case, it is advisable to retrieve the same contour from different images or to put more intermediary contours.From Landsat images collected in 1985 and 1986 we know that the water stage was higher than during the current decade but the method described here cannot be applied due to changes in dune and shoreline formation, which are inconsistent with the measure of ICESat during the last decade.These constraints and the need for close inspection may limit the applicability of this method, but in the case of Lake Poopó-with no gauging station and no radar altimetry coverage-even a partial bathymetry proves to be a major advancement.As expected, elevation precision seems to decrease (Table 2) with the number of profiles used in the calculations.It was not possible to compare absolute values of our hypsometry with previous works.Relative depth between surfaces of 500 km 2 and 2,250 km 2 in [3] seems to match the value of 2.25 m found in our study.Obtained volume and real water time series of Lake Poopó will greatly contribute to our understanding and knowledge of the changes occurring in the TPDS system.The outstanding vertical precision obtained in this study is mainly due to the total absence of vegetation close to the lake.This absence is due to water and terrain salinity and intense evaporation.In this study, only MODIS surfaces between 2000 and 2012 were used to retrieve continuous water storage time series, but there are 30 years of Landsat images in the archives that document past levels and storage variations.We believe that our method can be successfully applied in all similar regions with scarce vegetation, such as in the Middle East or for ephemeral lakes like Salar de Uyuni or Salar de Coipasa.As higher stages of the bathymetry can be easily approximated with existing DEM's, future studies with different methods and other remote sensing instruments will be concentrated on the lower part of the bathymetry of Lake Poopó where the last 25% of surface variations occur.The bathymetry and further modifications will be available through the LEGOS HYDROWEB service [49].

Figure 1 .
Figure 1.Lake Poopó, (a) geographical overview, (b) Ice, Cloud, and land Elevation Satellite (ICESat) ascending and descending tracks (yellow lines) over Lake Poopó.White circles represent zones where ICESat elevation is intersected with water contours.Maximum (red) and minimal (blue) water extents of our bathymetry are shown.

Figure 2 .
Figure 2. Typical ICESat profiles from northern Lake Poopó.Threshold black line 5 cm above water surface helps to exclude all points taken over water.(a) Example of profile taken during dry period (b) This profile shows many saturated points on the lake slope, which suggest crossing over surface with different optical properties (i.e., water).

Figure 3 .
Figure 3. ICESat profile is superposed on a Landsat image (white, red and green crosses).In addition, 2D plots show ICESat elevation values at corresponding latitudes.Water level is approximated by orange horizontal threshold line.Points above water (red) and saturated points below water (green) are shown.

Table 1 .
List of ICESat profiles with no visible saturation problems.

Table 2 .
Landsat cloud free and inundation free images used to retrieve water contours (date in format: sat_yyyy_mm_dd).Water contour elevation, elevation standard deviation (SD) and number of ICESat profiles used to calculate each elevation are given.