Comparative analysis of SRTM, TanDEM-X and UAV-SfM DEMs to estimate lavaka (gully) volumes and mobilization rates in the Lake Alaotra region (Madagascar)

Over the past decades, developments in remote sensing have resulted in an ever growing availability of topographic information on a global scale. A recent development is TanDEM-X, an interferometric SAR mission of the Deutsche Zentrum für Luftund Raumfahrt providing near-global coverage and high resolution DEMs. Moreover, ongoing developments in unmanned aerial vehicle (UAV) technology has enabled acquisitions of topographic information at a sub-meter resolution. Although UAV products are generally preferred for volume assessments of geomorphic features, their acquisition remains 5 time-consuming and is spatially constrained. However, some applications in geomorphology, such as the estimation of regional or national erosion quantities of specific landforms, require data over large areas. TanDEM-X data can be applied at such scales, but this raises the question of how much accuracy is lost because of the lower spatial resolution. Here, we evaluated the performance of the 12 m TanDEM-X DEM to i) estimate gully volumes, ii) establish an area-volume (A-V) relationship, and iii) determine mobilization rates, through comparison with a high resolution (0.2 m) UAV-SfM DEM and lower resolution (30 10 m) SRTM DEM. We did this for six study areas in the Lake Alaotra region (central Madagascar) where lavaka (gullies) are omnipresent and lavaka surface area changes over the period 1949-2010s are available for 699 lavaka. SRTM derived lavaka volume estimates were systematically too low, indicating that the SRTM DEM is too coarse to accurately estimate volumes of geomorphic features at the lavaka-scale (100 100 000 m). Lavaka volumes obtained from TanDEM-X were similar to UAV-SfM volumes for the largest features, whereas the volumes of smaller features were generally underestimated. To deal 15 with this bias we introduce a breakpoint analysis to eliminate volume reconstructions that suffer from processing errors as evidenced by significant fractions of negative volumes. This elimination allowed the establishment of an area-volume relationship for the TanDEM-X data that is within the 95% confidence interval of the UAV-SfM A-V relationship. Our calibrated area-volume relationship enabled us to obtain large-scale lavaka mobilization rates ranging between 18± 6 and 289± 125 ton ha-1 yr-1 with an average of 102± 41 ton ha-1 yr-1. These results indicate that current lavaka mobilization rates are two orders 20 1 https://doi.org/10.5194/esurf-2021-64 Preprint. Discussion started: 24 August 2021 c © Author(s) 2021. CC BY 4.0 License.


25
Over the past decades more and more advanced technology has become available for the assessment of surface topography: SfM (structure-from-motion) algorithms applied to UAV (unmanned aerial vehicle) imagery now allow centimeter-scale resolution, thereby revolutionizing the way we study earth-surface processes (Passalacqua et al., 2015;Tarolli, 2014;Clapuyt et al., 2016).
Obtaining high resolution DEMs from UAV-SfM still requires substantial fieldwork and is spatially limited due to the nature of the technology (Bangen et al., 2014). On the other hand, TanDEM-X is a remote sensing product with global coverage at 30 12 m resolution and, while being less detailed than high resolution DEMs, is a major step forward in comparison to the 30 m DEMs with a global coverage (Mudd, 2020).
This raises the question to which extent TanDEM-X imagery can be used to map three-dimensional morphological features requiring a higher degree of topographical detail over relatively large areas (> 10 km 2 ). One example is the use of remotely sensed data to map the process of gully formation that is known to significantly contribute to surface erosion (e.g. Poesen et al.,35 2003; Vanmaercke et al., 2021). Gully erosion mapping and monitoring was conventionally based on time consuming and spatially limited field surveys (Castillo et al., 2012;Evans and Lindsay, 2010;Guzzetti et al., 2012). More recently, however, high resolution DEMs have enabled the development of (semi-)automated gully-delineation and volume determination methods (Niculit ,ȃ et al., 2020;Evans and Lindsay, 2010;Perroy et al., 2010;Eustace et al., 2009;Liu et al., 2016), where TanDEM-X was already shown capable of automatically detecting gullies (Orti et al., 2019). 40 Not only the extent to which TanDEM-X data can be used to estimate gully volumes, but also its capability in establishing accurate area-volume (A-V) relationships is important to evaluate. This latter question is important since high resolution surface imagery from a multitude of sources and moments in time is now globally available. This imagery can be used to identify geomorphic features and estimate their surface area. A-V relationships then enable to obtain estimates of volume-changes over time when historical imagery from which areas can be derived is available. Work on landslide erosion has shown that 45 the establishment of area-volume relationships enables us to estimate sediment mobilization rates (i.e. the average annual volume of hillslope material displaced per unit area) over large spatial and temporal scales (e.g. Malamud et al., 2004;Larsen et al., 2010;Hovius et al., 1997;Guzzetti et al., 2012Guzzetti et al., , 2009. Furthermore, work on gully headcut retreat rates and associated sediment mobilization has indicated the importance of long measurement periods: large year-to year variations result in very large (>100%) uncertainties over short (<5 years) measuring periods (Vanmaercke et al., 2016;Frankl et al., 2013). 50 Here, we evaluate the performance of TanDEM-X to estimate gully volumes and to establish area-volume relationships by comparing estimates obtained from TanDEM-X with those obtained from a high resolution UAV-SfM and from a lower resolution SRTM DEM. We used the lavaka of the central highlands of Madagascar as a case study. Lavaka are amphitheater shaped 2 https://doi.org/10.5194/esurf-2021-64 Preprint. Discussion started: 24 August 2021 c Author(s) 2021. CC BY 4.0 License. gullies, with small outlets and are on average 60 m long, 30 m wide and 15 m deep (Cox et al., 2010;Wells and Andriamihaja, 1993). They are omnipresent in the central highlands, leaving the landscape filled with 'holes' ('lavaka' in Malagasy). Unlike 55 conventional gullies they typically lack surface feeder channels and tend to form on mid-slopes, broadening uphill trough headward erosion (Wells et al., 1991;Wells and Andriamihaja, 1993). While Madagascar is often claimed to experience amongst the highest global erosion rates due to the presence of lavaka (Milliman and Farnsworth, 2013;Randrianarijaona, 1983), the amount of sediment that is directly produced by lavaka is currently unknown.
The objectives of our study are therefore to evaluate the performance of TanDEM-X to i) estimate lavaka volumes, ii) 60 establish accurate area-volume relationships and iii) obtain a first estimate of lavaka mobilization rates. We derived lavaka volumes and mobilization rates for an existing dataset containing 699 digitized lavaka in six study areas in the Lake Alaotra region at three moments in time : 1949, 1969 and the 2010s. In a first step lavaka volumes were calculated from the DEM as the difference between a reconstructed pre-erosion surface and the current topography. Next, a lavaka area-volume relationship was established between the current lavaka areas and calculated volumes. Finally, this relationship was applied to the historical 65 dataset with lavaka areas in 1949, 1969 and the 2010s. This enabled to calculate lavaka volumes at each of these timesteps and the consequent derivation of volumetric growth rates and lavaka mobilization rates for each of the six study areas. This procedure was followed for a high resolution UAV-SfM DEM (0.2 m), the mid-resolution TanDEM-X DEM (12 m) and the low resolution SRTM DEM (30 m).

Study area and lavaka dataset
Six study areas (SA) of ca. 10 km 2 were selected in the northeastern part of the central highlands in the area surrounding Lake Alaotra (Fig. 1a). The lake is located in the seismically active NE-SW oriented Alaotra-Ankay graben structure and is surrounded by convex-shaped, deeply weathered (> 10 m) regolith-covered hillslopes that developed on the Precambrian crystalline basement (Kusky et al., 2010;Riquier and Segalen, 1949;Bourgeat, 1972;Mietton et al., 2018). The climate is 75 characterized by a distinct dry and wet season with the regular occurrence of tropical cyclones. The mature rounded hillslopes are covered with open grasslands and contain one of the highest lavaka densities of the country, with up to 14 lavaka km -2 (Cox et al., 2010;Voarintsoa et al., 2012). The lowland areas bordering the lake consist of swamps in the SE and vast rice fields elsewhere, producing the majority of rice for the country (Lammers et al., 2015). For the six selected study areas digitized lavaka polygons are available from orthorectified and georeferenced historical aerial images from 1949 and 1969 and 80 from recent (2011-2018 refered to as 2010s) satellite imagery (Maxar-Vivid-WVO2) (Fig. 1a, Table A1, Brosens (2020)). All shapefiles have WGS84-UTM 39S (EPSG: 32739) coordinates. The dataset contains the changes in surface area of 699 lavaka over the period 1949-2010s for SA1-5 and 1969-2010s for SA6. Each study area contains 50 to 173 lavaka, resulting in lavaka densities between 4 and 17 lavaka km -2 (Table A1). Lake Alaotra catchment shown on TanDEM-X DEM with hillshade (Krieger et al., 2007). The UAV-SfM DEM is available for study area 5 and 6 (red) (data collected June 2018). (b)-(d) Examples of the SRTM (Farr et al., 2007), TanDEM-X (Krieger et al., 2007) and UAV-SfM DEM (data collected June 2018) in SA6 located at 48°15'18.6"E 17°58'51.7"S with hillshade. Grey outlines indicate the digitized lavaka.

85
Lavaka volumes were determined from three digital elevation models with a range of horizontal resolutions. For two study areas a high resolution (0.2 m) UAV-SfM DEM was obtained from a field campaign in 2018. For all study areas the mid-(12 m) and coarse (30 m) resolution TanDEM-X and SRTM DEMs are available. All DEMs were transformed to WGS84-UTM39S (EPSG: 32739) coordinates.

90
For study area 5 and 6 ( Fig. 1a in red) a UAV-field survey was carried out in June 2018 to obtain a high resolution DEM.
In order to cover a large area during a limited amount of time with a high spatial resolution the post-processing kinematic 4 https://doi.org/10.5194/esurf-2021-64 Preprint. Discussion started: 24 August 2021 c Author(s) 2021. CC BY 4.0 License.
(PPK) georeferencing approach as developed by Zhang et al. (2019) was used. The UAV-images were directly georeferenced by using a RTK (real-time kinematic) receiver on the UAV which was connected to a RTK base station. This results in a robust and accurate alternative for georeferencing based on ground control points (GCP) (Zhang et al., 2019). Given that optical 95 acquisitions were georeferenced using RTK-GPS data, this high resolution surface can be considered as the reference of the 'true' elevation (Grohmann, 2018). In this study we therefore consider the UAV-SfM DEM as the ground-truth reference.

100
The RTK-receiver was connected to a single-board computer in order to synchronize the GPS time with the geotagging of the images. The RTK base station (Emlid Ltd) provided the positioning correction input. It was mounted on a tripod and located at a fixed position at the center of each study area during the flights. Flights were carried out at 200 m height at a speed of 8 m s -1 . Pictures were taken every 3.8 seconds resulting in an average ground sampling distance of 0.17 m ( Fig. 1e-f).
The raw RTK-GPS data from the receiver were corrected with the data from the base station and post-processed using the 105 RTKLib package (Takasu and Yasuda, 2009). A fix-and-hold method with 20°satellite elevation mask was used to correct the positions and geotag the images. The geotagged images were processed using Pix4D software using the default settings with the vertical and horizontal accuracy set at 0.5 m. For the generation of the DEM no surface smoothing or filtering was applied.
The resulting DEM for both study areas has a resolution of ca. 0.2 m (Fig. 1b).

110
The TanDEM-X (TerraSAR-X add-on for digital elevation measurements) mission was launched in 2010 by the public-private partnership between the German Aerospace Center (DLR) and EADS Astrium GmbH (Krieger et al., 2007). Its configuration consists of two synthetic aperture radar (SAR) satellites flying in close formation, thereby forming a large X-band single-pass interferometer. The resulting global DEM has a horizontal resolution of 0.4-arcsecond (ca. 12 m) and aims at 2 m relative height accuracy (Krieger et al., 2007(Krieger et al., , 2013. The final TanDEM-X DEM was published in 2016 and consists of data collected 115 between December 2010 and early 2015, where all land surfaces were imaged at least twice and up to 7 or 8 times in difficult terrain (Rizzoli et al., 2017, Fig. 1c).

SRTM (30 m)
The Shuttle Radar Topographic Mission (SRTM) resulted from a collaboration between the NASA, the National Geospatial-Intelligence Agency, and the German and Italian space agencies and collected images between 11 and 22 February 2000 aboard 120 of the Endeavour space shuttle (Farr et al., 2007). Data were collected from both X-and C-band radar interferometry, where the latter resulted in a first near-global DEM at 3-arcsecond resolution (ca. 90 m) that was released in 2005 (Farr et al., 2007).
In 2015 an enhanced void-filled 1-arcsecond (ca. 30 m) dataset was released (NASA-JPL, 2013), which was used in this study (Fig. 1d). For the African continent a 90% absolute geolocation error of 11.9 m and a 90% absolute height error of 5.6 m have been reported (Rodríguez et al., 2006).

Lavaka volume determination
Individual lavaka volumes were determined from the difference between the current surface and a reconstructed pre-erosion surface. This was done by developing an automated workflow in PyQGIS written in QGIS version 3.8.1 with GRASS 7.6.1 (code and example dataset available at https://doi.org/10.5281/zenodo.5155317). The automated PyQGIS workflow consists of six steps which are explained in detail below. The input data required to run the procedure is discussed in step 0.

STEP 0: Input data
Three input files are required to run the automated volume-procedure: i) a shapefile containing the digitized lavaka outlines, ii) a shapefile containing a pre-erosion surface polygon for each lavaka and iii) a DEM raster file. A manual delineation-procedure was followed to obtain the pre-erosion surface where a horseshoe-shaped polygon was drawn around each individual lavaka on the hillslope parts that were unaffected by erosion (Fig. 2a). This approach was pre-135 ferred over an automated pre-erosion interpolation (e.g. Evans and Lindsay, 2010) or interpolation based on the lavaka outlines for two reasons. First, digitized lavaka outlines from aerial imagery are often located on DEM pixels that already have lower values. This is especially the case for the coarser resolution DEMs due to surface smoothing (e.g. Each point is assigned the elevation value from the corresponding DEM-pixel.

STEP 3: Interpolate pre-erosion surface
The pre-erosion surface is obtained by interpolating between the pre-erosion polygon points. Two interpolation methods were used: i) Triangulated Irregular Network (TIN) and ii) spline interpolation. TIN interpolation is based on a linear method where triangles are constructed from the nearest neighbour points, resulting in non-smooth surfaces (Bergonse 150 and Reis, 2015; QGIS, 2020) (Fig. 2b). For the spline interpolation a regularized spline with tension algorithm was applied using the default settings (tension = 40, no smoothing), which enables the generation of curved surfaces in areas without data (Mitášová and Mitáš, 1993;Bergonse and Reis, 2015;GRASS, 2003) (Fig. 2c).

STEP 4: Calculate elevation difference
The current DEM is subtracted from the interpolated pre-erosion surface. The result is a difference raster with positive 155 values indicating a current surface that is lower than the reconstructed pre-erosion surface. Negative values indicate that the current topography is higher than the reconstructed topography, which is physically impossible.

STEP 5: Elevation difference clipped to lavaka extent
The lavaka extent, which is given by the digitized lavaka polygon, is clipped from the elevation difference raster. In this way a raster with the elevation difference over the lavaka area is obtained ( Fig. 2b and c). If the lavaka is smaller than one 160 pixel (0.04 m 2 , 144 m 2 and 900 m 2 for the UAV-SfM, TanDEM-X and SRTM DEM, respectively) the resulting raster is empty and no volume can be calculated.

STEP 6: Export results
The unique values report of the lavaka elevation difference raster is exported. It contains the unique elevation values, their count and dimensions of the raster pixels. These results are used to calculate the volumes of each lavaka.

165
From the exported values report the lavaka volume was calculated as the product of the elevation difference with the pixel area. Both positive and negative elevation differences occur, where a positive difference indicates that the current surface is lower than the reconstructed pre-erosion surface. Negative values indicate that the current topography is higher than the reconstructed topography. In principle negative values can result from two types of error: i) errors in the estimation of present heights and ii) errors due to the interpolation of the pre-lavaka surface area. It is not possible to distinguish between both error 170 types but it can be assumed that the first type of error was less important for the high resolution DEM. Given the presence of both positive and negative elevation differences three types of volumes (V [m 3 ]) could be calculated for each lavaka: i) the positive volume (V pos) calculated by summing only the positive elevation differences, ii) the negative volume (V neg) calculated from the negative elevation differences, and iii) the total volume (V tot) as the difference between the positive and negative volume. We also calculated the percentage of negative volume (V neg%) which is the fraction of the negative volume 175 over the absolute sum of the negative and positive volumes. If not further specified the term 'volume' refers to the positive volume.
Upon comparing the results from the different DEMs we considered the volumes from the UAV-SfM DEM as the groundtruth reference. This is, however, not entirely correct as these volumes also required a reconstruction of the original topographic surface which is subject to (unknown) errors, a universal problem in all research using reconstructed topography as a reference 180 (Bergonse and Reis, 2015). While the exact error cannot be determined since the original surface is unknown, the general performance of each DEM and interpolation method in determining lavaka volumes can be evaluated based on the relative proportion of the negative volume (V neg%): it can be reasonably assumed that a relatively large negative volume is associated with a relatively larger error in total volume.
The results from the different DEMs were compared in two ways: i) pairwise comparison using only the lavaka for which 185 the volume was determined for each DEM enabling a one-to-one comparison, and ii) using all the data in order to establish the most robust relationships calibrated on a higher number of observations. The number of lavaka for which the volume was determined depends both on the availability of the UAV-SfM DEM (only in SA 5 and 6) and the resolution of the DEM, where no volume could be calculated for the smallest features using the coarser resolution DEMs.  Krieger et al. (2007)) and SRTM (30 m, bottom, Farr et al. (2007)). (a) The digitized lavaka outline (grey), manually determined pre-erosion polygon (horseshoe-shaped polygon surrounding the lavaka) and current DEM are the three required inputs for the automated volumeprocedure (STEP 0). Random points are created in the pre-erosion polygon to which the DEM-elevation values are attributed (STEP 1-2).
The pre-erosion surface is then reconstructed by interpolating between these points. Two interpolation methods are tested: TIN (b) and Spline (c) interpolation (STEP 3). The outer grey polygon indicates the edge of the interpolated area. For both interpolation methods the elevation difference between the interpolated pre-erosion surface and current DEM surface is then calculated, which is clipped to the lavaka extent (STEP 4-5).

Establishing area-volume relationships 190
Establishing a relationship between lavaka area and volume enables to estimate lavaka volumes when only surface area information is available. Area-volume (typically landslides) and length-volume (typically gullies) relationships obey a power-law relationship V = aA b , where the predicted volume V for a given area A depends on the scaling exponent a and intercept b (Larsen et al., 2010;Frankl et al., 2013). A linear relationship is typically fitted on the log-transformed data in order to obtain 8 https://doi.org/10.5194/esurf-2021-64 Preprint. Discussion started: 24 August 2021 c Author(s) 2021. CC BY 4.0 License. equally distributed residual errors, resulting in a more robust fit: log(V ) = a + b log(A) (e.g. Guzzetti et al., 2009;Crawford, 195 1991). We established the relationship between lavaka area and volume by fitting a linear least-squares regression through the log-transformed data (base 10 log). When back-transforming the coefficients of the fitted linear relationship to a power-function a systematic statistical bias enters. This is accounted for by adding a bias-correction factor which depends on the variance σ 2 (Ferguson, 1986;Crawford, 1991): V = exp(a + 2.65σ 2 )A b . This correction assumes that the the residual errors of the fitted linear relationship are normally distributed with a mean of zero and variance σ 2 . The normal distribution of the residual errors 200 was tested using a Shapiro Wilk test (Shapiro and Wilk, 1965).

Lavaka volume growth and mobilization rates
where i indicates the most recent and j the oldest observation moment.
Lavaka mobilization rates (LM R [ton ha -1 yr -1 ]) give the amount of sediment that has been mobilized over a given period and area and were calculated for each study area. To obtain LM R, lavaka volumes were converted to mass using a dry 210 bulk density (ρ) of 1.5 ton m -3 (average dry bulk density from 2 m deep grassland soil corings in the Lake Alaotra region (Razanamahandry et al., submitted)). The sum of the lavaka masses of each study area was then divided by the length of the observation period and the surface area of the study area (A [ha]) : with N the number of lavaka in each study area.

215
Uncertainties on the calculated V GR and LM R were estimated by taking into account the uncertainties on the fitted a and b coefficients of the applied A-V relationship. This was done by running 10 000 Monte Carlo simulations with different a and b values. These values were randomly drawn from the normal distribution of their mean and standard error of the mean. Gaussian copulas were used to take into account the dependence of both coefficients (Frees and Valdez, 1998). From all 10 000 runs the mean and standard deviation were calculated.

TIN vs. spline interpolation
For the reconstruction of the pre-erosion surface both a spline and TIN interpolation were applied. The results of both interpolation methods were evaluated by comparing the resulting negative volume fractions (V neg% again non-significant for UAV-SfM (p = 0.07) and SRTM (p = 0.42). However, for TanDEM-X spline interpolation results in a significantly lower V neg% compared to TIN interpolation (p = 0.02, Fig. 3, Table A2). The interquartile range of V neg% is the lowest for the high resolution UAV-SfM DEM and increases with increasing DEM resolution, spanning over 90% for SRTM for both interpolation methods (3). The median V neg% is similar for the UAV-SfM and TanDEM-X DEM and is considerably higher for the SRTM DEM (Table A2).

235
While the difference between both methods is not significant in most cases, the systematically lower medians and visual interpretation of the resulting interpolated surface (e.g. Fig. 2) suggest that spline interpolation is the better option to reconstruct plan-curved surfaces which are omnipresent in our study area. This was also concluded by Bergonse and Reis (2015) who compared both techniques using a validation-based approach. They showed that spline methods result in smaller errors compared to linear (TIN) methods. Spline methods were shown to be better adjusted to a gully geomorphic context as they 240 allow curved surfaces in no data-areas, which is not the case for the linear interpolation of TIN (Bergonse and Reis, 2015). We therefore used the lavaka volumes derived from the spline-interpolated pre-erosion surfaces in the remainder of the analysis and manuscript.

Lavaka volumes
A direct pairwise comparison of the lavaka volumes obtained from the high resolution UAV-SfM DEM (0.2 m) and the coarser 245 resolution TanDEM-X (12 m) and SRTM (30 m) DEMs indicates that SRTM generally results in a large volume underestimation (Fig. 4a). Zero-volume instances represent lavaka for which the calculated total volume was negative and thus set to zero. TanDEM-X, on the other hand, results in similar volume-estimates as obtained by the UAV-SfM DEM, especially for the larger lavaka (> ca. 10 4 m 3 , Fig. 4a). While the absolute volume difference between UAV-SfM, TanDEM-X and SRTM increases with increasing lavaka volume (Fig. 4b), a different picture emerges when looking at the relative volume differences.

250
The relative difference is largest for the smallest lavaka and decreases when lavaka become bigger (Fig. 4c). Both absolute and relative volume differences are largest for SRTM, with strong volume over-and underestimations for the smallest lavaka and relative underestimations remaining above 80% for the larger features. For TanDEM-X large relative differences also occur for the smallest lavaka, however, they decrease to less than 20% for the largest features (Fig. 4c). The coarser DEM resolution of TanDEM-X and SRTM thus results in a systematic bias for all lavaka in the case of SRTM and for smaller lavaka in the case

255
of TanDEM-X. The less accurate and detailed topographic information of TanDEM-X and SRTM DEMs leads to a 'smoother' current topography and reconstructed pre-erosion surface, making smaller lavaka 'disappear'.
The smoothing effect of coarser resolution DEMs on landscape topographical representation is known to result in a reduced ability to capture more complex topography and geomorphic features (Thompson et al., 2001;Wechsler, 2007;Tarolli, 2014;Hengl, 2006). Furthermore, the optimal DEM grid resolution depends on the inherent properties and scale of the geomorphic 260 features under study (Tarolli, 2014;Hengl, 2006). Ideally a DEM should represent the properties in such a way that smallerscale features that will be filtered out do not harm the overall quality of the model outcome (Claessens et al., 2005). Our results clearly indicate that the 30 m resolution SRTM DEM is too coarse and filters out too many topographic details to accurately calculate volumes of erosional features at the lavaka scale (100 -100 000 m 2 ), where large errors remain even for the largest features (Fig. 4). Therefore the SRTM DEM will not be considered for further analysis.

Area-volume relationships
Area-volume relationships were established using the log-transformed data. The resulting pairwise relationships are fairly similar for the data obtained from the UAV-SfM and TanDEM-X DEM, although the smallest lavaka volumes are underestimated by TanDEM-X (Fig. 5a). This results in a lower intercept and corresponding higher scaling exponent for the TanDEM-X relationship, which is significantly different (outside of the 95% confidence interval) from the UAV-SfM relationship for the 270 smallest lavaka but within uncertainty for the larger features (Fig. 5a). A different picture emerges when considering the full dataset using all lavaka volumes for both DEMs: large volume underestimations for the smallest lavaka become apparent for volumes determined from the TanDEM-X DEM (Fig. 5b). These are caused by negative volumes that were calculated over (a fraction) of the lavaka area. For larger volumes, this discrepancy between both DEMs disappears, indicating that TanDEM-X is capable at accurately assessing lavaka volumes for features that exceed a given size. However, because of the large volume  While it is clear that the TanDEM-X DEM is too coarse to accurately predict lavaka volumes for the smallest features, this issue seems to disappear for the larger features. Therefore, we tried to identify the point below which the analysis based on 280 TanDEM-X suffers from errors in volume reconstruction as evidenced by negative volume pixels. This point was identified by determining the breakpoint in the V pos-V tot relationship when applying a broken-stick regression (Toms and Lesperance, 2003). The optimal breakpoint of the fitted piecewise linear spline model was selected based on the minimum sum of residuals of the fitted relationships using the SLM toolbox (D'Errico, 2021). This breakpoint is located at a positive volume of ca. 8000 m 3 and surface area of ca. 1900 m 2 (log(V pos) = 3.90, Fig. 6a). In a next step, we established a new area-volume relationship 285 for the TanDEM-X data containing only lavaka volumes larger than this identified breakpoint. This results in a close match with the fitted relationship based on the UAV-SfM data, where both regressions are within the 95% confidence intervals (Fig.   6b, Eq. (3) and (4)). Back-transforming the fitted linear log-transformed A-V relationships results in the following power-law still suffer from other resolution/smoothing effects. However, given the good correspondence between the results obtained with the high resolution UAV-SfM DEM and the results obtained with the TanDEM-X DEM it can be safely assumed that such errors will be mainly related the reconstruction of the original surface rather than to the lower resolution of the TanDEM-X DEM. Our V pos-V tot breakpoint method allows to exclude lavaka that suffer from evident volume reconstruction errors in an objective way. This method can furthermore be applied to regions where no high resolution DEMs are available for comparison, as the 300 breakpoint can be determined from the difference between V tot and V pos using the coarser resolution DEM. In our case, we can compare both resulting A-V relationships, confirming that this volume threshold results in a TanDEM-X A-V relationship that is within uncertainty of the high resolution UAV-SfM relationship (Fig. 6b). Furthermore, this threshold corresponds to the point where TanDEM-X volumes no longer deviate from volumes obtained from the UAV-SfM DEM (Fig. 4a). This seems to indicate that the largest errors in the volume estimates from TanDEM-X are contained in the percentage negative volume.

Lavaka volumetric growth and mobilization rates: 1949-2010s
By applying the established A-V relationships to the historical (1949 for SA1-5 and 1969 for SA6) and current (2010s) lavaka 315 areas (Table A1), volumetric growth rates (VGR) could be estimated (Eq. (1)). When using the UAV-SfM relationship (Eq. (3)) a mean and median growth rate of 907 ± 340 m 3 yr -1 and 265 ± 75 m 3 yr -1 are obtained, respectively. When applying the TanDEM-X relationship (Eq. (4)) these values are ca. 15% lower: 766 ± 117 and 228 ± 25 m 3 yr -1 for the mean and median, respectively. This deviation of 15% is, however, still within uncertainty of the estimates from the UAV-SfM DEM which have an uncertainty range of 37% resulting from the larger uncertainties on the fitted coefficients. While the scaling coefficient a is 320 identical for the TanDEM-X and UAV-SfM relationship, the intercept b for TanDEM-X is slightly lower (0.53 vs. 0.62, Eq. (3) and (4)). This indicates that small variations in the established coefficients can lead to relatively large differences in estimated volumetric growth and mobilization rates.
Volumetric growth rates (VGR) can be converted to lavaka mobilization rates (LMR) when the bulk density and size of the study areas are taken into account (Eq. (2)). Lavaka mobilization rates as derived from the UAV-SfM relationships range 325 between 18 ± 6 and 289 ± 125 ton ha -1 yr -1 in our six study areas (Table 1). LMR are again estimated to be ca. 15% lower when applying the TanDEM-X relationship (Table 1), resulting in LMR between 16±2 and 258±43 ton ha -1 yr -1 , within uncertainty of the UAV-SfM estimates.
While the estimates obtained using the UAV-SfM DEM and the TaNDEM-X DEM cannot be statistically distinguished, the lower estimates of VGR and LMR we obtained when using the TanDEM-X DEM are not unexpected. The underestimation of 330 erosion rates when coarser resolution DEMs are used was also reported by Claessens et al. (2005), who found that the highest landslide erosion and deposition was estimated for the highest resolution DEM, and estimated erosion rates systematically decreased when reducing the DEM resolution. This effect was attributed to the more detailed landscape representation for higher resolution DEMs. This likely also explains why our TanDEM-X based estimates are somewhat lower: the higher resolution UAV-SfM DEM will be able to capture more topographic details compared to the TanDEM-X DEM, resulting in slightly higher 335 calculated volumes and corresponding lavaka mobilization rates.   Lavaka mobilization rates (LMR) are positively correlated with the mean surface area of lavaka in the study area (spearman correlation coefficient r = 0.94, p = 0.02) (Fig. 7b). This can be explained by the positive correlation between lavaka area and volumetric growth rate (r = 0.27, p = 1e-10, Fig. 7a): larger lavaka mobilize more material. LMR also increase with increasing lavaka density (r = 0.94, p = 0.02, Fig. 7c), which is logical but is also partially explained by the positive correlation 340 between lavaka density and mean lavaka surface area (r = 0.89, p = 0.03, Fig. 7d). The main variations in LMR between our six study areas thus seem to depend mainly on the lavaka density and area distribution. While lavaka presence has been linked to seismic activity (Cox et al., 2010) and is typically associated with slopes ranging between 25 to 30°in the Lake Alaotra region (Voarintsoa et al., 2012), differences in surface growth rates could only be poorly linked to other environmental factors: the combined effects of the percentage bare surface area, distance to stream and distance to drainage divide could only explain In order to further evaluate the possible impact of fitting the TanDEM-X relationship on the larger features only (> 1900 m 2 ), we quantified the share of the total mobilized sediment that is provided by lavaka smaller than 1900 m 2 . From the relative cumulative sediment mobilization curves it is apparent that larger lavaka contribute most of the mobilized sediment (Fig. B2, note that the areas are plotted on a log-scale). Lavaka that are smaller than the identified threshold contribute 1.1% of the total 350 mobilized sediment in the study areas with the largest lavaka and up to 21.6% in the regions with smaller lavaka (Fig. B2).
This indicates that the share of smaller lavaka to the total amount of sediment that is mobilized is generally low in our study areas, therefore reducing the risk of erroneous estimates in the case where these smaller lavaka could not be used to establish the TanDEM-X based A-V relationships.
Our calculated sediment mobilization rates from direct lavaka growth observations over the period 1949-2010s (ca. 70 years) 355 range between 18 ± 6 and 289 ± 125 ton ha -1 yr -1 for six 10 km 2 study areas in the Lake Alaotra region with an overall average of 102 ± 41 ton ha -1 yr -1 . Only limited data is available that can be used to compare these estimates with. A sedimentation rate of 20 ton ha -1 yr -1 was obtained by Mietton et al. (2005) for the dammed Bevava lake which is located in the southeast of the Lake Alaotra catchment over the period 1987-2005. Lake Bevava has a catchment area of 58 km 2 with a lavaka density of 8 lavaka km -2 . We obtained a similar, though somewhat higher, erosion rate of 53 ± 19 ton ha -1 yr -1 for SA3, which has a similar 360 lavaka density of 9 lavaka km -2 (Fig. 7c, Table A1 and 1). It should be noted that our highest reported LMRs of 289 and 143 ton ha -1 yr -1 correspond to areas characterized by large lavaka and high lavaka densities (13 and 17 lavaka km -2 , Table A1, Fig. 7b). These lavaka densities are higher than the reported average of 6 lavaka km -2 for the southern part of the Lake Alaotra catchment (Voarintsoa et al., 2012). We therefore argue that these highest values should be perceived as maximum rates, where the rates of 18-53 ton ha -1 yr -1 obtained for regions with lower lavaka densities (SA 3, 5 and 6) will be more representative for 365 the wider Lake Alaotra region. Next to these recent short-term erosion rates, long-term catchment wide erosion rates obtained from 10 Be measurements have been reported for the central Malagasy highlands. These values range from 0.16 to 0.54 ton ha -1 yr -1 with the highest rates for the catchments with higher lavaka densities (max. 6 lavaka km -2 , Cox et al., 2009). Comparing these long-term erosion rates that integrate over tens of thousands of years with the recent short-term rates for the past 70 years indicates that current mobilization rates exceed the long term rates by two orders of magnitude. While not all mobilized lavaka 370 sediment will end up in the rivers, this still suggests that lavaka erosion has increased over recent time periods in the Lake Alaotra region.

Conclusions
Lavaka volumes were estimated as the difference between the current and interpolated pre-erosion surface for three DEMs with different spatial resolutions: SRTM (30 m), TanDEM-X (12 m) and UAV-SfM (0.2 m). Volumes estimated from TanDEM-X 375 are similar to those obtained from the UAV-SfM DEM for the larger features. Using the SRTM DEM results in strong volume underestimations, even for the largest features. This indicates that the SRTM DEM is not suitable to estimate erosion volumes for geomorphic features at the lavaka-scale (100 -100 000 m 2 ). TanDEM-X can be used for volume estimations, but shows a tendency to underestimate volumes for small lavaka, caused by a smoothing effect where complex topography and smaller geomorphic features cannot be accurately captured by its coarser resolution. An area-volume relationship, necessary for large 380 scale and past volume assessments, can be established using TanDEM-X or UAV-SfM data. However, developing a robust relationship based on the TanDEM-X data requires that observations for which the relative reconstruction error is too large are eliminated from the dataset. Here, we proposed and tested a method to identify a cut-off point below which volume estimations are clearly affected by processing errors as evidenced by negative volume estimates. This breakpoint is located at a lavaka volume of ca. 8000 m 3 , corresponding to an area of ca. 1900 m 2 . The proposed objective filtering to eliminate erroneous 385 volumes in the TanDEM-X estimates resulted in deviations in lavaka growth rates and mobilization rates that area ca. 15% lower compared to the respective UAV-SfM estimates and fall within the uncertainty boundaries of the latter. Our results thus indicate that the TanDEM-X DEM can be used to establish accurate area-volume relationships for erosional features at the lavaka-scale. As the proposed method does not depend on direct comparison with higher resolution DEMs, it can be applied to regions where only TanDEM-X is available. Over the period 1949-2010s a mean and median lavaka volumetric growth rate 390 of 907 ± 340 and 265 ± 75 m 3 yr -1 and lavaka mobilization rates varying between 18 ± 6 and 289 ± 125 ton ha -1 yr -1 were obtained. While our highest lavaka mobilization rates are likely limited to the most lavaka dense regions, our lower estimates are consistent with dam siltation rates, placing the current average lavaka erosion rate for the Lake Alaotra region at 20-50 ton ha -1 yr -1 . These rates are furthermore two orders of magnitude higher than earlier reported long-term erosion rates, suggesting a recent increase in lavaka erosion intensity in the Lake Alaotra region.

395
Code and data availability. TanDEM Figure B1. Example of near absence original surface topography. Example from study area 1 illustrating the near absence of the original surface topography (especially in the western part of the area) due to the dense presence of lavaka (grey outlines). Elevations from the TanDEM-X DEM with hillshade (Krieger et al., 2007).