Quantifying sediment mass redistribution from joint time-lapse gravimetry and photogrammetry surveys

. The accurate quantification of sediment mass redistribution is central to the study of surface processes, yet it remains a challenging task. Here we test a new combination of terrestrial gravity and drone photogrammetry methods to quantify sediment mass redistribution over a 1-km 2 area. Gravity and photogrammetry are complementary methods. Indeed, gravity changes are sensitive to mass changes and to their location. Thus, by using photogrammetry data to constrain this location, the sediment mass can be properly estimated from the gravity data. We carried out 3 joint gravity-photogrammetry surveys, once 20 a year in 2015, 2016 and 2017 over a 1-km 2 area in southern Taiwan featuring both a wide meander of the Laonong River and a slow landslide. We first removed the gravity changes from non-sediment effects, such as tides, groundwater, surface displacements and air pressure variations. Then, we inverted the density of the sediment, with an attempt to distinguish the density of the landslide from the density of the river sediments. We eventually estimate an average loss of 3.7 ± 0.4×10 9 kg of sediment from 2015 to 2017, mostly due to the slow landslide. Although the gravity devices used in this study are expensive 25 and need week-long surveys, new instrumentation progresses shall enable dense and continuous measurements at lower cost, making this method relevant to improve the estimation of erosion, sediment transfer and deposition in landscapes


Introduction
The reliable quantification of sediment mass redistribution is critical to the understanding of surface processes (Dadson et al., 2003;Hovius et al., 2011;Morera et al., 2017) and has significant implications for studies in tectonics (Molnar et al., 2007;Steer et al., 2014;Willett, 1999), climate (Peizhen et al., 2001;Steer et al., 2012), human activities (Horton et al., 2017;Torres et al., 2008;Hare et al., 2008).When studying underground processes, especially groundwater, it is common to simplify the redistribution to a one-dimension problem.The groundwater level variations Δh are the main observable and gravity effects are computed using a Bouguer plate (2  Δℎ).This simplification is necessary because it is usually impossible to monitor lateral groundwater redistribution, and the assumption of little lateral change remains appropriate for homogeneous aquifers.The groundwater level variation can be assumed constant over the entire aquifer.Such an assumption is not valid for surface processes because sediment builds complex three-dimensional bodies.But sediment mass redistributions occur at the ground's surface; thus, they are accessible to accurate location methods such as photogrammetry (Eltner et al., 2016;Niethammer et al., 2012;Schwab et al., 2008).Combining accurate geometries with gravity variations can thus enable proper mass estimations.Fig. 2 illustrates the use of time-variable gravimetry to quantify sediment mass redistribution at the earth's surface.In the simplest case, when considering each ground element as a point-mass, the total change of gravity Δg measured between t0 and t1 is: where Δgi is the vertical component of the gravitational change at each element i (i ranging from 1 to N = 28 in Fig. 2b) considered as a point-mass (Fig. 2c) of mass mi located at a distance ri from the gravimeter, and G is the universal gravitational constant.Note that the gravitational attraction of any element decreases with the square of the distance between this element and the site where gravity is measured, so that the distance of the mass redistribution can be a strong limiting factor to measuring significant gravity changes.Note also how the angle θ between the point mass and the site where gravity is measured contributes to the gravity effect.The gravity effect is maximum when the point mass is at the vertical of the site, negative if above the site, positive if below.If the point mass is exactly at the horizontal with the gravimeter sensor, then the gravity effect cancels.The effects of the angle and the distance are shown in Fig. 2d and 2e, for a general case and for one actual site of the survey, respectively.Point-mass simplification is ideal to grasp the concept of gravimetry but it is not suitable for precise quantifications that are the aim of this study.All gravity modelling will thus be done using rectangular prism modelling (Nagy, 1966), which is the most appropriate way to compute the gravity effect of surface changes measured by photogrammetry.
After introducing the study area, we describe the gravimetry and photogrammetry surveys that we conducted, together with our data processing workflow.We then show the results of both methods and interpret them jointly in order to retrieve the mass of sediment redistributed in this area from 2015 to 2017.We eventually discuss the benefits and limits of this method.

Study area
The joined gravimetry-photogrammetry survey was set in southern Taiwan, at the Paolai village, next to the Laonong River (Fig. 1).The gravity network contains one site, AG06, for absolute gravity (AG) measurement and nine sites, BA01 to BA09, for relative gravity (RG) measurements.During the 2017 survey, all sites but BA02 were located to cm accuracy using Global Navigation Satellite System (GNSS) enhanced by real-time kinematic (RTK) technique.The exact location of BA02 could not be measured due to the unexpected storage of concrete blocks, referred to as dolosse, placed on the river shore to protect it from erosion.This dolosse storage also covered BA03 and BA04 but those two sites could still be measured.The gravimetric effect of the dolosse was estimated and removed from the measurements.
The first reason for choosing this location is that time-lapse absolute gravity surveys have been done at AG06 since 2006, in the frame of the Absolute Gravity in the Taiwan Orogen (AGTO) project.This project made it possible to measure, for the first time, sediment mass redistribution using time-lapse absolute gravimetry and showed that significant sediment transfers occurred around Paolai (Kao et al., 2017;Mouyen et al., 2013).Indeed, this site experiences vigorous sediment transfer processes powered by heavy rains brought by tropical cyclones (typhoons) and monsoonal events, especially in May to August (Chen and Chen, 2003).The heavy rains destabilize the slopes of the Taiwanese high mountains, triggering landslides and debris flows (Chiang and Chang, 2011).This occurs on a regular basis: 5 to 6 typhoons make landfall in Taiwan every year (Tu et al., 2009), mostly between May and September.The most remarkable event was the 2009 Morakot typhoon, which produced the worst flooding in the last 60 years in Taiwan and up to 2777 mm of accumulated rainfall (Ge et al., 2010) and triggered 22 705 landslides with a total area of 274 km 2 (Lin et al., 2011).Landsliding, which can also be triggered by regional active tectonics, is the main process supplying sediment to rivers in Taiwan (Dadson et al., 2003;Hovius et al., 2000).
The second reason for choosing this location is practical.This location offers a stabilized path made of concrete on the southern bank of the Laonong River, where the relative gravity benchmarks could be properly set, on stable and sustainable sites, and easily accessed for measurements.Also, a continuous GNSS station, PAOL (latitude: 23.10862°, longitude: 120.70287°, elevation: 431 m), is co-located with the AG06 pillar and is maintained by the Institute of Earth Sciences-Academia Sinica (IES-AS, 2015).This makes it possible to precisely take into account gravity changes only due to ground vertical displacements.
In this area, both the Laonong River and the landslide (Fig. 1) are susceptible to sediment transfers.The gravimetryphotogrammetry survey is set up to focus on these processes.Note that what we call the river (plain black line contour in Fig. 1) is the active channel bed that includes emerged alluvium.During yearly measurements, the water extent of the river only covers a fraction of this area, even if the period 2015-2017 saw some higher water level and larger extents, especially during large floods.

Methods
This section introduces the two main methods used in this study: gravimetry and photogrammetry.Gravimetry is sensitive to masses and their distribution, while photogrammetry is here a geometric measure of the ground surface, hence of the sediment distribution.Therefore, combining gravimetry and photogrammetry removes the geometric ambiguity inherent to gravity measurements and allows to focus on sediment masses.This combination is done through a least-squares inversion to determine sediment density, that is a mass per unit of volume.

Time variable gravimetry
Gravity was measured at 10 sites (Fig. 1) in 2015, 2016 and 2017, always over two days in November, since the climatic conditions during this month are usually suitable for gravimetry fieldwork (e.g.no typhoon nor heavy rains, reasonable temperatures).By measuring gravity during the same period of each year, we also expect to minimize hydrological effects, which have an strong annual periodicity in this area (Chen and Chen, 2003).Absolute and relative gravity surveys were done in parallel, on the same days.
Relative gravity measurements were done using a Scintrex CG5 Autograv (serial number 167).The measurement principle is to assess length variations of a spring holding a proof mass between different times and places, using a capacitive displacement transducer, and convert them into gravity variations (Scintrex Ltd., 2010).The instrument is levelled at each site and repeats 90-seconds measurements continuously.We stop measurements when gravity readings repeat within 3 µGal (1 µGal = 10 - 8 m s -2 ), while the internal sensor temperature remains stable.This usually takes 10 to 15 measurements, that is 15 to 23 minutes, although up to 25 measurements were required in some rare cases.Only the latest measurements, when gravity readings are stable, are used in the gravity network adjustment, to estimate the drift.Indeed, relative gravimeter measurements are subjected to an instrumental drift, which is corrected using the software Gravnet (Hwang et al., 2002).Inferring this drift requires re-measuring one or more sites within a few hours.In this study, all surveys start and end at AG06, which is also revisited up to four times during the survey, together with other relative gravity sites (Appendix A).In addition, ambient temperature alters gravity measurements at a rate of -0.5 Gal °C-1 (Fores et al., 2017).This effect was taken into account before adjusting the instrumental drift of the gravimeter.
The absolute gravity measurements were done using a Micro-g FG5 (serial number 224), which monitors the drop of a freefalling corner-cube in a vacuum.During its free fall, the positions and times of the corner-cube are precisely assessed using laser interferometry and an atomic clock (Niebauer, 2015;Niebauer et al., 1995).One measurement takes ~12 hours and consists of 24 sets of 100 test mass drops started every 30 min (one drop every 10 s).Measurements are always done overnight, when anthropogenic seismic noise and temperature variations are lower than during day time.A laser problem in the FG5 prevented us from measuring absolute gravity in 2017.This is compromising since the measurements at BA01-BA09 need an absolute reference to be compared with the previous survey in 2016 of these sites.In order to interpret the 2017 relative gravity survey, we decide to estimate the AG06 absolute gravity value in 2017 as the mean of the values measured in 2014, 2015 and 2016, that is 978 713 845.1 ± 3 µGal (Fig. 3).We believe this is a reasonable approach because there were no major climate or tectonic events between November 2016 and November 2017 around the Paolai region.In addition, by repeating the gravity surveys at the same time of the year, in November, the gravity difference due to hydrological changes are at the few-microgal level, which is about the size of the errors in the relative gravity measurements.The hydrological conditions are described in more details in the next paragraphs.For the 2017 estimated absolute gravity value, we arbitrarily set the standard deviation to a value about three times larger than usual at this site.The AG06 values in 2016 and 2017 are quite similar, less than 0.5 µGal difference.Although absolute gravity measurements at AG06 started in 2006 and repeated once a year except from 2011 to 2013, it is not possible to use these older data for the estimation of the 2017 value.Indeed, in 2009, Typhoon Morakot and its subsequent massive landslides reset the whole area.The gravity offset between November 2008 and November 2009, i.e. before and after Morakot, is about 30 µGal and is due to large sediment redistribution in this area (Mouyen et al., 2013).
Sediment redistribution due to Morakot was such an exceptional event, with a significant impact on gravity, that it cannot be included in the extrapolation of the 2017 gravity value.The measurements from 2009-2010 were not used either because too much reconstruction work was ongoing at that time, taking out debris from the river, thus interfering with natural sediment redistribution.
To focus on sediment mass redistribution, other sources responsible for gravity changes must be removed from the gravity time series.Here, these effects are the tides, air pressure variations, polar motions, vertical ground motions and hydrology.
These corrections are detailed in the next paragraph and summarized in Table 1.
Solid earth tides are computed using TSOFT (Van Camp and Vauterin, 2005) using tidal parameters from Dehant et al. (1999), referred to as WDD.Ocean tide loading effects are computed using the FES2004 model (Letellier et al., 2004) with the Ocean Tide Loading provider (Bos and Scherneck, 2003).Polar motion effects are computed using the International Earth Rotation and Reference System Services parameters and the Absolute Observations Data Processing Standards (Boedecker, 1991).
Atmospheric effects, that is gravity changes due to air masses, are corrected using local barometric records done at a continuous weather station located ~12 km west of Paolai (station C0V250) and an admittance factor of -0.3 µGal hPa -1 (Merriam, 1992).Solid Earth tides, ocean tide loading and atmosphere loading are corrected before the drift adjustment of the relative gravity measurements, because they can have significant effects over a few hours, that is while the relative gravity survey is done.Not correcting them would bias the drift estimation by mixing gravity changes due to the instrumental drift with those due to tides and atmosphere.Vertical displacements of the ground also change the gravity, because the gravity measured at any place on the Earth's surface depends on the inverse of the square of the distance between this site and the Earth's center of mass.
Therefore, if the site is uplifting (further from center of mass) or subsiding (closer to the center of mass), it will have a lower or higher gravity value, respectively.This effect is corrected using continuous GNSS time series recorded at AG06 (the GNSS site PAOL is co-located with AG06, Fig. 4) and assuming a theoretical ratio Δg/Δz = -2 µGal/cm (Van Camp et al., 2011), where Δg is the gravity change and Δz is the elevation change, at the same location.Between 2015 and 2017, the ground uplift at AG06 is about 1.3 cm yr -1 .That is a large uplift rate, explained by the active mountain building processes at work in Taiwan, where up to 1.9 cm.yr -1 uplift is measured (Ching et al., 2011) .Although the relative gravity sites are between 300 and 500 m from the PAOL permanent GNSS, we apply the same uplift correction to these sites as to AG06.Indeed, tectonic uplift is a regional feature and can be assumed constant over a few hundred meters, unless an active fault or more cross the area.But there is no evidence for such a fault in Paolai.
We also correct the effect of hydrology, which deforms the earth surface at the global scale and changes the groundwater mass attraction at local scales, near the gravimeter.This correction relies on global hydrological models.We consider two of them in this study: 1. the Global Land Data Assimilation System Version 2 (GLDAS-2) forcing the Noah land surface model (Rodell et al., 2004) and 2. the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2, Gelaro et al., 2017).
The gravitational effect due to each of these models is provided by the EOST loading service (Boy, 2015;Petrov and Boy, 2004).Unlike the other corrections, the hydrological correction may suffer large uncertainties because of (1) the complexity of hydrological processes, (2) the difficulty to measure groundwater and (3) its large effect on gravity (Jacob et al., 2009;Longuevergne et al., 2009;Pfeffer et al., 2013).Indeed, the effect of GLDAS-2 and MERRA-2 on gravity predict up to 20 µGal of seasonal amplitude in the hydrological signal, with sometimes large differences, up to 10 µGal, between the different models (Fig. 5).Nevertheless, surveying in November appears to be a valuable way to decrease the hydrological impact on the gravity data, since this effect is lower than 3 µGal, when considering either of the two hydrological models.We use the average hydrological effect from GLADS-2 and MERRA-2.We arbitrarily set an uncertainty of 5 µGal to this correction (Table 1), to account for possible bias in the models.
We also correct the effect of the dolosse set in 2017, which is only significant at BA03 and BA04.These structures, located above BA03 and BA04, were responsible for an artificial decrease of gravity of about 15 µGal.Their gravitational effect is computed using the dolosse's geometry measured by the photogrammetry and rectangular prism method (computation details in Appendix B).Given the uncertainty of this correction process, we add an arbitrary 5 µGal uncertainty on the gravity changes measured at BA03 and BA04 during the 2017 survey.
The last non-sediment effect on gravity is the actual Laonong River: water with density ρw = 10 3 kg m -3 .The photogrammetry measures the river surface but not its depth.Also, at each survey, the river did not follow the same path across the active river bed.Therefore, neither the volume of the river nor its variation can be estimated, which prevents us from a gravity correction for the river.

Photogrammetry
The purpose of the photogrammetry survey is to generate high resolution digital surface models (DSM) in 2015, 2016 and 2017, at the moment of the gravity surveys, to quantify the sediment volumes changes between each survey.

Equipment
The photogrammetry survey was done with an unmanned aerial vehicle (UAV), commonly known as a drone, which is commonly used in morphotectonics studies (Chang et al., 2018;Deffontaines et al., 2017Deffontaines et al., , 2018)).The UAV is a modified already-available Skywalker X8 fixed-wing aircraft reinforced by carbon fiber rods and Kevlar fiber sheets (Figs.6a and 6b).
To generate a high-resolution DSM, orthorectified mosaic images, and a true 3D model, the UAV is equipped with a Sony ILCE-QX1 global shutter camera and a 16 mm SEL16F2.8 lens.

Survey design and execution
The UAV is launched by hand, then flies, takes photos, and lands autonomously by using a pre-programmed flight plan.The autopilot system is composed and modified from the open source APM (Ardupilot Mega 2.6 autopilot) firmware and open source software Mission Planner (Oborne, 2010), transmitted by ground-air XBee radio telemetry.
The UAV was flown between 300 and 500 m above the ground level.It covered an area of about 15-20 km 2 with about 15-20 cm ground sampling distance (GSD) in one single flight mission (Table 2).Repeated adjacent photographs were kept for at least 85 % endlap and 50 % sidelap.The GPS location of the acquired image is recorded on the auto-pilot log.Each UAV flight missions took about 90 min.On average 13 ground control points (GCP) per survey and 11 check points (CKP) were measured in the field to control and verify the quality of the datasets (Fig. 6c).The GCP are pre-existing or painted benchmarks on the ground.The CKP are permanent features such as bricks with patterns, road sign and zebra crossing.The GCP and CKP coordinates were measured by Virtual Base Station (VBS) RTK GPS and RTK GPS at least 3 times each.The mean vertical error and the root mean square are 1.5 cm and 4.2 cm, respectively.

Photogrammetry processing and results
The images acquired by the UAV, their positions and the GCP positions are processed using Pix4Dmapper (Pix4D, 2020) in order to generate a DSM, with a grid spacing of 50 cm based on multi-view stereo, which aims at reconstructing a complete 3D object model from a collection of images taken from known camera viewpoints (Furukawa and Ponce, 2010;Seitz et al., 2006).Images and camera parameters such as focal length, principal point, radial distortion, tangential distortion, aspect ratio and skew are auto-adjusted for each image during the processing.The point cloud densification and DSM parameters set specifically for Pix4Dmapper are summarized in Table 3.
The quality of the DSM is calculated by comparing its elevation with those of the CKP at the same coordinates.The error of the dataset is denoted as Table 4, where the data shows that the mean of the error compared with field survey about 0.11-0.20 m and with standard deviation of 0.334 to 0.622 m.

Survey results
The results of the gravity and photogrammetric surveys are summarized in Figs.7 and Table 5.The largest gravity changes occurred between 2015 and 2016, with most sites showing an increase of more than 30 µGal.In contrast, the gravity decreased at most sites from 2016 to 2017.When measured above the redistributed masses, increase and decrease of gravity correspond to gain and loss of masses, respectively.Qualitatively, this is in agreement with the corresponding changes observed in the digital surface models (DSM) in the active bed channel showing higher sediment thicknesses, thus a gain of mass, from 2015 to 2016 and large surfaces of lower sediment thicknesses from 2016 to 2017.Over the time period 2015-2016, the top of the landslide is actively eroded, up to 46 m, while its toe displays significant sedimentation, up to 33 m.The active river bed shows a mixed-pattern of erosion and sedimentation, between -1.19 and 1.21 m on average, possibly resulting from the migration of the river braids.In contrast, over the time period 2016-2017, the landslide displays mostly erosion, up to 39 m, while the river bed continues to display a mixed-pattern of erosion and sedimentation, between -1.17 and 1.08 m on average.
The gravimetric and photogrammetric techniques show large changes in gravity and topography, which demonstrate active processes of sediment mass redistribution in the river and on the slow landslide.In the next section, we combine these results to assess the mass of sediment redistributed from 2015 to 2017.Note that we focus the DSM analysis on the area bounded by the black line in Figs.7d and 7e, which is restricted to the landsliding zone and the river.

Joint analyses of the gravity and photogrammetry data
Using the DSM, we build rectangular prisms with horizontal sides of 0.5 m, i.e. at the scale of the DSM resolution, and a vertical side as high as the elevation at the time of the corresponding surveys, i.e. bottom at 0 m and top at the surface elevation.
Among the three (2015, 2016 and 2107) photogrammetric surveys, the 2017 survey has the smallest surface extent.Its limits are thus used to cut the 2015 and 2016 photogrammetric surveyed areas, so that all DSMs cover the exact same area.The total mass of redistributed sediment equals the change in volume between each survey multiplied by the density of the sediment.
We use the gravity data to assess this density using an inverse modelling approach.Note that since gravity decreases with the square of the distance between the measurement site and the mass location, we can bound our analysis to the area covered by the photogrammetric surveys without biasing the analysis.Indeed, using the wider 2015 and 2016 survey coverages, we find that extending our working area in the north-south and east-west directions by steps of 100 m does not alter the gravity changes computed at each sites by more than 1 %.
We design three inversion cases to retrieve the densities of the redistributed materials, using a least-square criterion.These cases are independent from each other and aim at increasing the amount of possibly different densities for comparison.Thus we invert:  Case 1: The average density ρ of the material redistributed during the surveys.
 Case 2: The density of the sediment in the river ρr and the density ρl of the material in the landslide.
 Case 3: The density of the sediment in the river ρr 1615 from 2015 to 2016 and ρr 1716 from 2016 to 2017 and the density ρl of the material in the landslide.
Here we will solve an over-determined problem, where we have more observations (20 gravity differences over the three years) than variables to estimate (density, three at most, in case 3).However, gravity observations are too few and unevenly distributed over the study area to try to invert the density at each pixel (more than 4 million) of the photogrammetry survey.In practice, the matrix representation of this system is (e.g.Hwang et al., 2002)  +  = (2) where the design matrix A, vector of unknowns X and vector of observations L are defined as (5) and V is the vector of residuals (X and V are to be determined by the least-squares method).In matrices A and L, dg is the gravity variation that is modelled (dgmod) or observed (dgobs) between 2016 and 2015 (1615) or between 2017 and 2016 (1716) at every site (BA01… AG06).The modelled gravity change can be computed for the material in the river (dgmod,r) or in the landsliding zone (dgmod,l).This matrix representation is given for the inversion case 3 and can be simplified for cases 1 and 2.
The design matrix A is built thanks to the photogrammetry surveys, from which we identify the river and the landslides as well as their respective volume changes.Therefore, knowing also the position of the gravity sites, we compute each element of A using a gravity modelling by rectangular prisms methods (Nagy, 1966) and an arbitrary density equal to 1.The actual density can be inverted by where A T is the transpose of A. The weight matrix P is diagonal, and its elements are the inverse of the gravity uncertainties at each site i.The residuals  =  −  are used to compute a posteriori variance of the unit weight where n is the number of gravity observations and u the number of unknown densities.The uncertainties of the inverted densities are the square root of the diagonal element of the a posteriori covariance matrix of X   =  0 2 (  ) −1 (8) The inverted densities for each case are summarized in Table 6.Cases 1 and 2 return similar densities.Case 3 returns a noticeable difference between the densities of the sediment in the active bed channel for the 2015/2016 or 2016/2017 surveys.
A first hypothesis for this difference could be that the composition of the redistributed sediment has changed over the years of the study, for instance because material comes from landslides that occurred in terrain with different densities.A second hypothesis is that the water content of the sediment varies, changing the effective density of the sediment as measured by the gravimeters.We do not have enough data to favour one of these hypotheses but we will discuss the possible influence of water on our density estimates in section 6.1.Uncertainties on the landsliding materials densities (case 2 and 3) are higher than those of the river sediment, likely because they are further from the gravity sites than the river sediment.As seen in equation 1, the further the redistributed masses are, the lower are their gravitational effects.
We described in section 3.1 the impossibility to remove the effect of the Laonong River (water) from gravity changes due only to sediment redistribution.As a workaround, here we simply assume a constant river depth of 1 m, which corresponds to rough field estimates.Then, we map the surface limits of the river from the optical images taken by the UAV during each survey.
The height h of the river surface is given by the photogrammetry results.The river is then divided into prisms covering the river area, with sides of 0.5 m, upper face at elevation h and lower face at elevation h-1, since the river is 1-m deep.We then compute the gravity effect of the river at each site of the network (except BA02, which position is unknown).This effect is removed from the gravity observation and the average density inversion (case 1) is run, giving  = 1.7 ± 0.1 × 10 3 kg m -3 and RMS = 9.7 µGal.This represents a decrease of 11% relative to the density  = 1.9 ± 0.1 × 10 3 kg m -3 given in Table 1 and also relative to the mass budget in Fig. 10.These values are to be taken with caution since we do not know the exact geometry of the river, its depth in particular.In the rest of the study, we will only work with the density given in Table 1.
For comparison, during the 2017 survey, we evaluate the in situ densities of the materials in the active river bed and at the bottom of the landslide, at 22 sites (Fig. 8), also using photogrammetry (Appendix C).Estimating in situ density is timeconsuming and demanding.It is done here only for comparison purposes; it is not required for the inversion.Indeed, joint gravity-photogrammetry estimates an average in situ density over the surveyed area.In contrast, in situ density measurements are done at discrete locations over an area made of heterogeneous material.The in situ densities range from 1.2 to 2.7×10 3 kg m -3 and are spatially heterogeneous, illustrating the variety of materials carried by the river and the landslide.Despite the limited and spatially uneven sampling points, we obtain an average density (2.0×10 3 kg m -3 ) consistent with the spatiallyintegrated density inverted from the gravity and photogrammetry data (1.9×10 3 kg m -3 ).It is interesting to note that the density in the lobes of the main, fresh, landslide (Fig. 8: the samples the most in the north, around 219500 m on the eastern axis) are among the lowest density measured.This landslide material is sourced from rocks, mostly slates that may have a higher density, about 2.7×10 3 kg m -3 (Ho, 1986).The landslide broke them and stacked them into a rough, unorganized pile, less compact than the original material.As a result, using the average density of the rocks in this area would overestimate the mass.
The final comparison of the measured gravity and the computed gravity in cases 1, 2 and 3 is given in Fig. 9.The largest misfits are at BA05 and BA06 during the 2016-2017 period, for which gravity changes are underestimated by 19 and 15 µGal, respectively.Possible explanations for these two misfits are: a wrong site location, an error in the gravity data, an error in the DSM data or local but large hydrological effects, not accessible at the scale of the global hydrological models we used.
However, we could not narrow our search down a specific issue at BA05 and BA06.At the other sites, the pattern and amplitude of the gravity observations is rather well explained by the modelling.Note that in Fig. 9b, the gravity modelled at most sites seems to need a small offset of -3 µGals to fit within the uncertainty of the observations.This may show that our absolute gravity estimate for the 2017 survey (Fig. 3) is wrong by 3 µGals.
Multiplying the inverted densities (Table 6) with the volume changes computed from the DSM changes, we can eventually compute the mass of sediment that was redistributed between two surveys, for each inversion cases (Fig. 10).Since the inverted densities are similar in each case (Table 6) and the volume changes estimated from photogrammetry are identical, the estimated masses (volumes times density) are also similar in each three case.The difference mostly lies within the uncertainty of these estimates.In our mass estimation, we also differentiate the top and the toe of the landslide, because the top of the landslide mostly experiences erosion, while its toe undergoes both erosion and sedimentation.This helps to unravel how the sedimentation and erosion processes are distributed over the slow landslide.
In the river only, we observe that the mass of sediment redistributed between each survey is similar.The river gained between 0.61 and 0.83×10 9 kg and lost between 0.58 and 0.74×10 9 kg.Thus, the mass loss is about 4% and 12% less than the mass gain, resulting in a quasi-balanced budget that is within the uncertainty of the mass estimations.The time variability of the sediment mass budget is dominated by the landslide, which causes larger mass redistributions (up to 4×10 9 kg) and loss-togain ratios.A significant mass loss occurred between 2016 and 2017, which is ~15 times larger than the mass gain.Between 2015 and 2016, both erosion and sedimentation are significant at about 2 to 3×10 9 kg, which are rather balanced.A likely hypothesis is that we mainly observe a transfer of sediment from the top of the landslide, where 2.1 ± 0.4×10 9 kg of material was eroded toward its toe, where 1.9 ± 0.4×10 9 kg accumulated (average mass from the three cases).Overall, from 2015 to 2017, the area has lost about 3.7 ± 0.4×10 9 kg of sediment.Note that this landslide occurs over several years, not in one quick event, probably as the consequence of the erosion by the meandering Laonong River.

Implications for sediment transfers in active landscapes
Our results highlight how landscapes react to landsliding and how they evolve after a large perturbation such as the 2009 Morakot typhoon.Between 2015 and 2016, the activity of the Paolai slow landslide mostly consists in transferring about 2×10 9 kg (about 1×10 6 m 3 ) of material from the landslide top to the landslide toe over roughly 100 to 200 m.After 2016, a significant event of erosion of the landslide occurs, with more than 3×10 9 kg of sediment removed, including most of the sediment previously deposited on the landslide toe.This corresponds to a particularly rapid evacuation of the sediment, especially in the alluvial context of the Laonong River, that is consistent with predictions obtained with a morphodynamic model by Croissant et al. (2017) for bedrock rivers.It is likely that the position of the landslide in the outer bank of a meander has favoured sediment export efficiency.Despite this landslide activity, it is quite remarkable that the Laonong River roughly maintains a neutral sediment budget over 2 years, between 2015 and 2017, in the vicinity of the landslide.This means that the river mainly acts as a sediment transfer zone and that river incision and sediment evacuation occurring along the river is balanced by the sediment delivery occurring by the supply of landslide materials.This sediment supply may originate from the several large landslides triggered in the Paolai area by the 2009 Morakot typhoon (Lin et al., 2011) and the following massive sediment aggradation along fluvial valleys up to 10, 30 and even possibly 100 m (DeLisle and Yanites, 2018;Hsieh and Capart, 2013).Our results would thus suggest that the Laonong River has not yet recovered from this aggradation phase and that the landscape is still perturbed by the aftermath of Morakot typhoon, even 8 years after its occurrence.This exceeds the relaxation time of 6 years observed after the 1999 Chi-Chi earthquake using river suspended load (Hovius et al., 2011) but is consistent with longer evacuation time-scales of coarser, non-suspended, materials (Yanites et al., 2010).Typhoon-triggered landslides occur every year in Taiwan, and global warming may intensify this process (Chiang and Chang, 2011).This could also build and maintain long-term sediment sources within the Taiwan range, which will keep supplying sediment into rivers even long after the Morakot-induced sources have been completely flushed.

Perspectives from recent progresses in gravimetry
In this study we take advantage of the intense surface processes occurring in Taiwan to jointly analyse both time-variable gravity and photogrammetry measurements.Indeed, the amplitude of the sediment mass redistribution guarantees that we were able to measure significant gravity changes and, most importantly, surface elevation changes.Nevertheless, for rivers experiencing large and dynamic sediment mass redistributions that remain hidden beneath the water level, photogrammetric data would not bring any constraint on the density inversion.One would thus only be able to rely on the gravity measurements, leading to a non-uniqueness problem, since both the density and the location of the redistributed sediment would have to be inverted.To better deal with this issue, we suggest two improvements to our gravity survey: 1. Set a denser network of gravity sites, ideally with a mesh structure.Indeed, more measurements, evenly distributed, mean more constraints on the density inversion.
2. Set this network closer to the mass changes to increase the gravity signal.The best option would be to locate the gravimeters directly beneath the river bottom.Fig. 11 shows that for such gravimeters, the average gravity change would increase from 31 to 50 µGal between 2016 and 2015 and from 13 to 61 µGal between 2017 and 2016.
This suggested survey design implies that gravimeters are set permanently over the time period of the project, as they would not be easily accessible.Such a setup of buried permanent gravimeters is presently impossible to realize with CG5 or any other contemporary relative or absolute gravimeter, but remains realistic at a few-years' timescale.Indeed, a new generation of relative gravimeters is rising from the use micro-electro-mechanical systems (MEMS) technology, characterized by a significantly smaller size and lower price (Liu and Pike, 2016;Middlemiss et al., 2016Middlemiss et al., , 2017)).These shoebox-sized devices could be used to set permanent and dense gravity networks in rivers.However, setting persistent gravity networks in areas experiencing vigorous sediment transport shall require deeper practical thinking.Gravity changes densely sampled over the river will make it possible to retrieve the sediment mass redistribution using gravity inversion methods (e.g.Camacho et al., 2011) further constrained by the geometry of the river and the depth of the relative gravimeter.In addition, as relative gravimeters suffer from instrumental drift, this buried permanent network should be run in parallel with either permanent absolute measurements, which has recently become possible thanks to quantum gravimeters achieving 1 µGal repeatability (Ménoret et al., 2018), or to slowly drifting superconducting gravimeters (Hinderer et al., 2015).Therefore, ongoing progress in the development of terrestrial gravimeters may open new opportunities for quantifying the mass of sediment redistributed by surface processes.Another interest for having such a permanent gravity network is to monitor the dynamics of the sediment mass redistribution at timescales shorter than one year, since the sediment concentration in Laonong River varies across the year (Fig. 12a).

Continuous sediment transport estimation
The method and perspectives introduced so far aim at quantifying the mass of sediment redistributed by an event with large sediment transport ability, such as a landslide or a high river discharge.The time step of this quantification depends on how long these events take to redistribute the sediment in a way that significantly alter the gravity measured at each site by, e.g., >10 µGal, as an indicative change.Nevertheless, the best solution is to set a permanent gravity network, so that any rapid sediment mass redistribution can be recorded.Fig. 12a shows that the largest sediment concentration recorded at LiuGui station is 5000 ppm (mass fraction), when the river level increased by 1.6 m.
For the hypothetical permanent buried gravity network (Fig. 11c), we compute the gravity effect of a river level change of 1.6 m, which covers the entire area of the active bed channel of the Laonong River.The 5000-ppm sediment concentration means there is 5 kg of sediment in 1000 kg of river fluid, hence 995 kg of pure water.In this framework, and assuming that the density of the sediment is 2×10 3 kg m -3 , we can compute the density change due to rising sediment concentrations, until 10 6 ppm, meaning the river is made of sediment only.The gravity variation solely due to a shift from 0 to 5000 ppm of suspended sediment is about 0.17 µGal on average over each site.This cannot be properly distinguished from the main gravity change due to the rising river water level.In fact, the suspended sediment concentration would need to increase by about 3×10 5 ppm to change the gravity by at least 10 µGal (Fig. 12b with the bedload set to 0 cm).This corresponds to a concentrated debris-flow, nearly 8 times more concentrated than the threshold of 4×10 4 ppm used for debris-flow definition (Dadson et al., 2005;Lin and Chen, 2013).However, sediment is also transported on the river bed, as bedload, and its variation must be added to the variation of suspended sediment concentration to estimate the effect of sediment discharge variations on gravity.We have no measurement of this bedload component for Laonong River but measurements in another catchment of Taiwan showed that 50% of the cumulative mass of the bedload was built by rocks which diameter is 15 cm (D50 = 15 cm) and D90 = 62.5 cm (Cook et al., 2013).Therefore, we compute the effect of adding homogenous bedload layers of up to 60 cm thickness and density 2×10 3 kg m -3 to suspended sediment variations (labelled curves in Fig. 12b).This addition generates about 50 µGal of gravity variation, which would be clearly identifiable in the gravity measurements.This computation gives an order of magnitude of the gravity change expected from time-varying suspended and bedload transport.It shows that continuous timevariable gravity could quantify changes of sediment discharge if the sediment concentration rises by at least 3×10 5 ppm without bedload, or if the bedload increases to a thickness of at least 12.5 cm, under the assumption that only gravity changes above 10 µGal are significant.More accurate predictions of gravity effects require knowing the proportionality relation, if any, between the suspended and bedload component, as well as local hydrogravity models.Again, this 10-µGal threshold is linked to the accuracy of todays' gravimeter but ongoing progresses and interest in time-variable gravimetry may fuel the development of devices with higher accuracies.

Conclusion
This study shows that the mass of sediment redistributed by rivers and landslides can be estimated by combining time-lapse gravimetric and photogrammetric measurements.Focusing on the Laonong River, in southern Taiwan, we estimate that about 3.7 ± 0.4×10 9 kg of sediment was removed from 2015 to 2017 around our study site.This sediment loss is mainly due to a slow landslide moving from one year to another.The sediment budget (i.e. the difference between sedimentation and erosion) within the river is close to zero, although more surveys should be carried out to identify longer-term deposition or erosion in this area.The average sediment density obtained with this method (1.9 ± 0.01×10 3 kg m -3 ) is similar to the average sediment density measured in situ across the flood plain (2.0×10 3 kg m -3 ).The new method introduced in this paper has the advantage of being able to directly sense the mass of sediment and can benefit a wealth of studies on surface processes, which require quantitative estimates of sediment mass redistribution.Although time-variable gravimetry remains a rather expensive method with demanding survey constraints, it has undergone promising progress in the recent years.One is the significant miniaturization of the devices, using inexpensive MEMS technology (Middlemiss et al., 2016), and the other is the realization of permanent absolute gravimeters, using cold-atom interferometry (Ménoret et al., 2018).Such new tools could be further used without photogrammetry, for rivers where most of the sediment transport is hidden under the water.If the suspended and bedload transport are significant enough, measuring the instantaneous sediment discharge could also become a reasonable 480 project.

Site selection:
We select 20 sites in the active river channel, which are accessible by walking.We try to find sites where materials are different, some of them being close to each other, to better grasp the variety of the material in the channel.Nevertheless, we also try to have a spatially even sampling.In a few places, two samplings are done at the same horizontal position but at the surface and then deeper.All site positions are recorded with a hand GPS (about 3 m accuracy).

Material sampling
At each start, we first distribute several benchmarked rules all around the place that will be sampled (Fig. C1a).Several pictures are taken to cover the sampling site and several benchmarked rules at a time.Pictures should overlap each other.We then dig a hole of about 30-40cm depth and radius, paying attention to not move any of the benchmarked rules.The excavated material M is put in a bucket of known mass and weighed using a hook-hanging weight machine.Then another set of pictures is taken to cover again the benchmarked rules and the hole just dug.The only difference between the pictures taken before and after is the hole.

Photogrammetry:
The benchmarked rules make a common reference between the pictures taken before and after the hole is dug.They are transformed into clouds of points in 3D (Fig. C1b), one representing the original surface, the other representing the dug surface.Thus, subtracting these two surfaces returns the surface of the hole, from which the volume V of the hole is computed (Fig. C1c).

Density computation:
The density at the sampling site is then M divide by V.

Figure 4 :
Figure 4: Ground vertical displacements time series at PAOL GNSS station, co-located with AG06, provided by the GPSLAB database (IES-AS, 2015).Solutions are computed in the IGS08 reference frame(Rebischung et al., 2012).The time of each joined 830

Figure 5 :
Figure 5: Hydrological effect on gravity at AG06, estimated from global hydrological models GLDAS2 and MERRA2.

Figure 6 :Figure 7 :Figure 8 :
Figure 6: a) UAV, modified Skywalker X8. b) Close-up on the central compartment of the UAV, where the camera and lens are mounted.c) Map of the ground control points and checkpoints (numbered according to Table 4) with the shaded topography in the background.The gravity sites are also shown for reference.

Figure 9 :
Figure 9: Comparison of the observed (blue) and modelled gravity changes for the densities inverted in cases 1, 2 and 3 (red, yellow and purple, respectively).Each case is slightly offset horizontally for legibility.No gravity is estimated at BA02 since its location is unknown.

Figure 10 :
Figure 10: Estimation of the mass of sediment redistributed between 2016 and 2015, and between 2017 and 2016 in cases 1, 2 and 3 (red, yellow and purple, respectively; same color code as in Fig. 9).The mass estimation is shown for the river, the landslide and 870

Figure 11 :
Figure 11: Gravity changes expected at new sites located 5 m beneath the river (red), compared with those measured at BA01-BA09 (blue), for the same sediment mass redistributions as a) between 2015 and 2016 and b) between 2016 and 2017.The new sites are in fact BA01 -BA09 translated 140 m in the north-east direction, as illustrated in c).

Figure 12 :
Figure 12: a) River level and sediment concentration of Laonong River, measured at LiuGui station, about 20 km downstream from Paolai.The highest sediment concentration (5000 ppm) is reached in summer 2017, when the river level increased by about 1.6 m.Data are freely available at the Taiwan WRA (Taiwan Water Resources Agency, 2019) .b) Estimated gravity changes at the buried network (Fig. 11c) as a function of the variations of the suspended sediment load and of increasing amounts of bedload-transported sediment.The bedload fraction is considered here as a homogenous layer of 0 to 60 cm thickness (labelled on each curve) and density 2×10 3 kg m -3 .The river becomes a debris-flow when its suspended sediment concentration goes beyond 4×10 4 ppm.The 5000 ppm level is shown as a reference.Note that the gravity sign is negative because the mass is increased above the gravimeters.

Figure A1 :
Figure A1: (a) Map view of the relative gravity network with the link between each site for the 2015 survey.(b) Gravity readings on the CG5 at each site as a function of time.(c) Histogram of the residuals after the drift adjustment.

Figure A3 :
Figure A3: (a) Map view of the relative gravity network with the link between each site for the 2017 survey.(b) and (c) Gravity readings on the CG5 at each site as a function of time.(d) Histogram of the residuals after the drift adjustment.

Figure
Figure C1: a) Picture of the hole taken with references scales and benchmarks.Several pictures were thus taken before and after the hole was dug.b) 3D cloud of the points mapping the surface of the hole.c) Computation of the volume bounded by the hole and the former surface of the ground, before the hole was done.920

Table 4 :
Results of field check points and the elevation difference (Δz) of photogrammetric DSM for each survey