Quantifying water storage change and land subsidence induced by reservoir impoundment 1 using GRACE , Landsat , and GPS data 2

12 The construction of hydropower dams is a common strategy to support a country’s increasing need for 13 electricity and river water management for industry and agriculture. Although the hydrological and 14 geophysical impacts of water relocation are usually assessed prior to impoundment, their accuracy is 15 generally limited due to the lack of in situ observations, especially in a remote area. This study 16 presents a workflow to quantify the terrestrial water storage change (∆TWS) and land subsidence 17 induced by a reservoir’s water impoundment using multiple satellite observations (GRACE, Landsat), 18 land surface models (CABLE, GLDAS, NCEP, ECMWF), and GPS data. The study site is the Bakun 19 Dam, located in Sarawak, Malaysia, which is the largest hydropower dam in Southeast Asia. 20 Commencing operation in late 2010, the dam induced a change of water mass and lake surface area 21 that was clearly observed by GRACE and Landsat observations, respectively. During the 17-month 22 impounding period (from August 2010 to December 2011), GRACE observed a dramatic increase of 23 approximately 200 mm equivalent water height, while Landsat detected an increased lake extent of 24 around 600 km. In this paper, a forward model is developed to determine the increased water surface 25

level corresponding to GRACE observations, estimated to be about 120 m. In contrast to GRACE, the 26 TWS derived from land surface models cannot capture the increased ΔTWS, due to the lack of 27 reservoir routing algorithms in the models. In addition, the land subsidence was calculated using the 28 disk load model constructed based on the GRACE-derived lake level and Landsat-derived lake extent; 29 the result is validated with the GPS data from BIN1 station, located at the western coast of Borneo. 30 The commencement stage of the Bakun Dam induces the large-scale land subsidence, which causes 31 the GPS-BIN1 station to subside by ~9 mm, and move toward the Bakun Lake by ~4 mm. 32

Introduction 42
Construction of fresh water reservoir is one of the potential solutions to secure the country 43 development under the risk of climate change and projected increases in population (Ehsani et al., 44 Lettenmaier & Famiglietti, 2006;Zhou et al., 2018). In contrast to LSM, GRACE observes the entire 79 water column, including all natural and anthropogenic processes, leading to a more complete and 80 accurate record of the estimated ∆TWS. Another advantage of GRACE is the sensitivity associated 81 with mass concentration. A gravity signal of a concentrated mass over a small area can be recovered 82 (Castellazzi et al. 2016), demonstrated by, e.g., Longuevergne et al. (2013), Castellazzi et al. (2018), 83 making GRACE a potential measurement for monitoring the storage variability of the local reservoir. 84 A significant water impoundment during the dam commencement generally induces a displacement of 85 the Earth's surface (Boy & Chao, 2002). The displacement can be computed using the GRACE SHC  (Han, 2016). 89 Due to the limited spatial resolution (or limited maximum harmonic degree, Nmax in the spectral 90 domain) of GRACE, the approach is only valid when the spatial extent of surface mass change is 91 sufficiently large (e.g., >250 km resolution). However, the GRACE results may be not adequate for 92 determining the mass change of a small spatial scale and the associated elastic surface deformation. 93 Such a phenomenon can be understood as a GRACE truncation error, which commonly leads to a 94 discrepancy between the GRACE and GPS's vertical displacement estimate (e.g., Khan  scenes over Bakun Lake, the yearly-averaged NDWI is used for the analysis. 157 The lake volume can be approximated using a lake filling model associated with the derived NDWI 158 and digital elevation model (DEM). Three DEMs are considered in the lake volume calculation, the 159 Shuttle Radar Topography Mission (SRTM, Rabus et al., 2003), the Advanced Spaceborne Thermal

GRACE data 165
The GRACE twin satellites measure variations of the Earth's gravity field every month using a 166 combination of several measurements onboard, e.g., K-band ranging, accelerometer, attitude, and 167 orbital data (Bettadpur, 2012). As the gravity variation at the monthly time scale is dominated by the 168 hydrological mass changes, the GRACE gravity data are often used to compute ∆TWS (Wahr et al., 169 1998)

GPS measurements 223
The GPS records at station BIN1 (Fig. 1a) are obtained from two data processing centers, the 224 series. To reduce the plate tectonic effect (e.g., Fu and Freymueller (2012); ), the associated long-term trend of up, north, and east component is firstly fitted (using Eq. (A1)) and 233 removed from its time series. Finally, the long-term mean is subtracted from the time series to obtain 234 the variation of the GPS position. 235 236

Methods 237
The processing workflow is shown in Fig. 2. The satellite images obtained from Landsat 7 and 8 are 238 used to delineate the Bakun Lake's extent as well as to approximate the lake's level. GRACE data and 239 LSMs are used to estimate the ΔTWS, and the comparison between GRACE and LSMs results is 240 conducted to evaluate a unique benefit of GRACE data in the controlled environment. A forward 241 modeling of the GRACE-derived ΔTWS anomaly during the water impoundment period is also 242 performed to approximate the water level. The lake's outline obtained from Landsat and the lake's 243 level estimated from the GRACE forward modeling are then used to construct the disk load model. 244 The land subsidence is assessed from the disk load model and the GRACE SHC estimate, and the 245 comparison between both results together with the GPS records will confirm the effectiveness of the 246 disk load model application over the Bakun Lake. 247 The CSI is computed as: 259 (2) 260 where is the Landsat pixel's index, is the number of total pixels, and 261 simulated pixels fail to detect water otherwise . 264 The CSI values vary between 0 and 1 where 1 indicate that the simulated water level produces an 265 identical flood extent as the observation. 266 The accuracy of derived NDWIs is assessed using a referenced image from the Earth Science Data

Derivation of lake's water level using TWS forward model 295
The forward model ( ) associated with the GRACE-derived ∆TWS can be used to determine a water 296 level of the Bakun Lake. The lake basin function ( ) is firstly designed by setting a value of a 297 uniform water level ( ) in the lake while setting outside value as zero. Then, the Gaussian smoothing 298 ( ) is applied to the basin function ( = ( )) to simulate the water storage at the GRACE spatial 299 resolution. The variance ratio, = ∑ ( − ) 2 2 ⁄ ( is the considered grid cell, and is the 300 number of grid cells), is used to find the optimal water level parameter from the GRACE 301 observation ( ). The spatial resolution of GRACE is obtained from the correlation length, estimated 302 from the empirical covariance function ( ) (see Eq. (3 -5) of Tscherning and Rapp, 1974), where 303 is the spherical distance. The correlation length is defined as the distance ( ) where the ( = 0) 304 decreases by half. 305 It should be noted that the GRACE-derived (equivalent) water level variation includes all 306 hydrological components, e.g., surface water, soil moisture, groundwater, and the estimate might be 307 different from the measured reservoir level (e.g., from gauge data) that reflects only the surface water 308 component. 309 310

Development of disk load model 311
The disk load is designed as a cylinder with a given radius and thickness. A number of discrete loads 312 are placed along the outline of the water body (e.g., delimited from the satellite imagery data), and the 313 total displacement is computed by summing the contributed responses from all disk loads. Following 314 this concept, the surface mass load ( , kg/m 2 ) as a function of the angular distance ( , radian) away 315 from the load center can be formulated based on a series of Legendre polynomials as follows: 316 where is the angular disk radius (radian), is the Legendre polynomial of degree (unitless), and 320 is the disc thickness (meter), which is defined as a water depth in this study. The vertical ( ) and 321 horizontal ( ) displacements (meter) caused by the disk load are formulated as:

Inundated area and lake volume from Landsat data 330
The yearly average NDWIs derived from the Landsat data between 2007 and 2016 are shown in Fig.  331 3. After a careful inspection, we found that the NDWI values greater than 0.2 represent the surface 332 water, while the values <0.1 generally include clouds and vegetation. Therefore, we define a surface 333 water pixel for NDWI > 0.2 as a safe threshold. The derived NDWI has an overall accuracy of ~95 % 334 evaluating against the ESDR-GFCC data (Sect. 3.1). We found that the accuracy is reduced when the 335 NDWI threshold is smaller than 0.2, which is likely caused by the contamination of cloud and 336 vegetation pixels in the calculation. From Dam began its operation (Fig. 3e), and the lake extent has no significant change from 2012 onward 340 ( Fig. 3fj). The averaged NDWI between 2012 and 2016 is then used to delineate the lake boundary 341 (polygon). The inundated area can be calculated by integrating the area of all water pixels inside the 342 lake polygon, as shown in Fig. 3k. The inundated area is substantially increased by 12 times (from The Landsat data and DEM are used to approximate lake volume. The water pixels derived from the 346 averaged NDWI between 2012 and 2016 is defined as the observation (Fig. 4a). Figs. 4bd 347 demonstrate the simulated inundated area after filling the lake area with the averaged water levels of 348 approximately 46 m, 60 m, and 73 m, respectively. With increasing water level, the flood firstly 349 emerges around the downstream of the Balui River and progresses upstream through Rumah Baka and 350 Rumah Kulit, respectively. The maximum CSI value is obtained when the simulated water level of ~ 351 75 m (corresponding to 48 Gt of total water mass) is used (Fig. 5), and it results in the best agreement 352 with the observed flood extension (Fig. 4a). Comparing with SRTM, the average water level is 353 reduced by ~ 2m, and increased by ~1 m when AW3D30 and ASTER-GDEM2 are used, respectively. 354 The difference is ~ 1.5 m or 1.8 % of the approximate water level. 355 feature, and the correlation length of ~275 km is found in all solutions (Fig. 7b). The forward model 391 shows a similar spatial feature to the observation around the dam location (Fig. 7c), and its intensity is 392 proportional to the parameter . Considering only the averaged ∆TWS inside the Bakun Lake, the 393 minimum variance ratio ( ) is found when ≈ 120 m (Fig. 7d). We assume that the error of GRACE-394 derived ∆TWS follows a normal distribution (described by zero mean and standard deviation of 52 395 mm), and the uncertainty of the approximate water level is estimated based on an error propagation 396 method. This corresponds to the water level uncertainty of ~27 m. The TWS forward modeling 397 approach is applied to all monthly GRACE-derived ∆TWS to approximate the water level variations of the Bakun Lake, and the water level time series is displayed in Fig. 2 (associated with the right 399 axis). 400  During the reservoir impoundment, the impounded water is first distributed through the soil and 435 groundwater zones, filling pore space until the maximum saturated TWS capacity is reached, and the 436 exceeded capacity results in surface water. The required filling capacity (before the appearance of the 437 lake) can be approximated by CABLE, whose land cover parameters are publicly accessible. We 438 found that the CABLE-estimated TWS between January 2010 and July 2010 (period prior to 439 impoundment) is stable at ~10.40 ± 0.35 m, which is lower than the saturated capacity by ~0.83 m. 440 Note that the obtained value is at the model resolution. We compute the required filling capacity at the 441 Bakun Lake by applying the TWS forward model in a similar manner as GRACE and obtained the 442 water level of ~ 34 ± 7 m that needs to be filled as a foundation of the Bakun Lake. 443 The summation between the Landsat-derived lake level (75 ± 1.5 m) and the required filling capacity 444 (34 ± 7 m) shows a reasonable agreement with the water level approximated by GRACE (109 ± 7 m 445 Vs. 119 ± 27 m). 446 The solid Earth deforms elastically in response to surface water load (Farrell, 1972). An elastic Earth 449 model is taken to compute the three-dimensional surface deformation induced by the Bakun Lake load 450 using the GRACE (truncated) SHC data as well as the disk load model of the reservoir. The disk 451 model is essentially determined with Landsat inundated area and constrained by GRACE-derived 452 ∆TWS for water level changes, as described previously. The Bakun Lake load is spatially localized 453 demanding spherical harmonic expansion to the degree considerably higher than what GRACE data 454

provide. 455
The Bakun Lake disk model is composed of multiple disks with variable radii covering the flooded 456 area, as shown on the inset map in Fig. 9a. The total area of the disk loads is ~600 km 2 , corresponding 457 to the area derived from the Landsat data (Sect. 3.1). The averaged monthly water level derived from 458 GRACE is used to constrain monthly changes of uniform thickness of the Lake disks (Sect. 3.3). The 459 vertical displacement is calculated from the disk model using Eq. (9) with max = 40,000, and the 460 result is shown in Fig. 9a. It is seen that significant subsidence (~200 mm) occurs beneath the lake a 461 result of the elastic deformation with the lake load and the amount of subsidence rapidly decreases 462 with distance from the lake. The subsidence is computed to be around ~9 mm at GPS site BIN1. The 463 vertical displacement is also computed using the GRACE SHC data using Eq. (4). The difference in 464 the GRACE-derived vertical displacement between the pre-dam and post-dam period is shown in Fig.  465 9b. The subsidence derived from the GRACE SHC data is markedly smaller than the disk model (e.g., 466 ~3 mm vs. ~200 mm around the lake) due to truncation of GRACE SHC data. The load deformation 467 signal spreads (leaks) over a wider area around the lake resulting significantly smaller subsidence at 468 the BIN1 site. Such leakage effect is considerable if the source load is spatially confined. 469 distributions is very similar (Fig. 11). Only a slight difference is seen near the upstream of the Bakun 491 Lake (Fig. 11a, b). The deformation profile between the BIN1 station (point A) and the Bakun Dam 492 (point B) shows a clear difference near the Bakun Dam, but hardly noticeable beyond 10 km away 493 from the lake (Fig. 11c). 494 the BIN1 station (point A in Fig. 11) and the Bakun Dam (point B in Fig. 11) are shown in Fig. 12. 498 Similar to Fig. 10c, the difference is noticeable near the Bakun Lake, where PREMsoft model are converged after ~40 km away from the lake, and the effect of using different Earth's model is 501 trivial at BIN1 station. Wang et al. (2012) and Wahr et al. (2013), also report similar behavior. 502 [ Fig. 12] 503 The computed vertical displacements from the disk load model and the GRACE SHC data are 504 compared with the GPS measurements at the BIN1 station after the long-term mean is removed (Fig.  505   13). The total subsidence is measured by two GPS solutions; GRACE, and the disk model, as shown 506 in Table 4. The total change is calculated by multiplying the estimated linear trend between August 507 2010 and December 2011 (see the displayed slopes in Fig. 13) with the impoundment period (e.g., 17 508 months). The GPS data present that the BIN1 station indicates subsidence of approximately 10 -11 509 mm (or about 7 -8 mm/year) during the impoundment period. The GPS solutions show a significant 510 agreement, with correlation value (ρ) of 0.9 (Fig. 14). The disk model is consistent with the GPS time 511 series (~ 9 mm subsidence with ρ = 0.77). Due to SHC truncation, the GRACE-derived vertical 512 displacement shows 6 -7 times smaller deformation at BIN1, and the ρ value is approximately 13 % 513 lower than the disk model. 514 [ Fig. 13] 515 [ Table 4] 516 [Fig. 14] 517 The horizontal displacements computed from the disk load model and the GRACE SHC data are 518 compared with the GPS measurements (Fig. 15). The disk load model predicts the deformation 519 southward by ~2.4 mm and eastward by ~3.3 mm (see Table 4), or 4.1 mm toward the Bakun Lake 520 during the impoundment period (Fig. 16). The motion of the GPS-ULR solution shows a reasonable 521 agreement with the disk load model or 3.8 mm toward the lake. The GPS-UNR solution only agrees 522 with the disk load model (and the GPS-ULR solution) in the north component, while it moves away 523 from the lake in the east-west direction (Fig. 15c). It is unclear why the east component of the GPS 524 time series is less consistent between two solutions. The inconsistent use of reference frame, SHC data fails to capture the elastic deformation by the lake load due to the omission error. Instead, 527 the computed deformation from GRACE SHC data is likely the deformation caused by the seasonal 528 migration of degree-1 load (geocenter motion), as also seen in Han (2016). The GPS horizontal 529 measurements also present the coherent seasonal motions synchronized with the GRACE results. Borneo island by ~200 mm in equivalent water height. This is 2-3 times larger than the seasonal water 543 storage change predicted by the LSMs, which do not include the Bakun lake storage but are 544 constrained principally by soil moisture and groundwater storages. We find the corresponding water 545 level of ~120 m increasing during the impoundment, which is the composite of surface water (70 %) 546 and underground water components (30 %). In total, GRACE observes ~71 Gt of fresh water that has 547 been stored in the Bakun Lake since 2011. This is comparable to nearly twice the annual ice mass loss 548 in Alaska glacier and more than the total annual loss of the western Antarctica ice sheet (Bamber, 549

2018). 550
Our developed disk load model of the Bakun Lake constrained by the GRACE and Landsat data 551 predicts land subsidence about 200 mm at the lake and ~9 mm around the nearby coastal area. The magnitude of the subsidence is around twice the global mean sea level rise trend (Chen et al., 2017). 553 As such, the relative sea level rise near the west coastal Borneo is supposed to be tentatively 2-3 times 554 faster during the impoundment period, and the impoundment effect should be considered in the 555 estimation of the sea level variation in the coastal zones (Nicholls and Cazenave, 2010). Such load 556 displacement estimate is in a good agreement with the actual GPS records at the coastal town Bintulu 557 (BIN1). The reservoir impoundment causes the GPS site to move toward the lake by ~4 mm, observed 558 from the disk load model and GPS-ULR solution. The magnitude of the horizontal displacement is 559 approximately two to three times smaller than the vertical component, which is in line with the 560 finding of Wahr et al. (2013). We emphasize here that the subsidence estimate of the Bakun Lake is 561 the result of the elastic response. Subsidence process could link to an inelastic response of Earth 562 material, e.g., a variation of aquifer system (Amelung et al. 1999;Erban et al., 2014), and our 563 approach is unlikely suitable for such an application. 564 This study presents a framework for merging the strength of multiple remote sensing data and 565 numerical modeling sources. The advantage and limitation of different data learned from this study 566 are summarized in Table 5. GRACE possesses a clear strength of observing the total mass variation 567 that cannot be captured by other observations. However, the direct use of GRACE data cannot provide 568 a complete gravity signal due to the truncation of SHC. The effect is severe for a load associated with 569 short wavelength, which can cause systematic misfits with the validated data (e.g., Fig. 11, see also Fu 570 et al., 2013). A numerical model can be developed to utilize the strength of GRACE and Landsat data 571 in an iterative manner and is used to recover the full spectral band of the gravity signal. As a result, 572 GRACE is capable of detecting a variability of localized mass distribution over the area as small as 573 600 km 2 . The size of this area is considered the smallest for which GRACE has been utilized for 574 hydrological application, e.g., ~10 times smaller than Lake Nasser (Longuevergne et al., 2013), and 575 Lake Urmia (Tourian et al., 2015). However, the mass variability of the Bakun Lake is as large as (or 576 larger) than Lake Nasser and Lake Urmia, which explains the sensitivity of GRACE over the Bakun

Conclusion 582
We developed and demonstrated a technique to quantify the total water storage change and land 583 subsidence caused by the impoundment of the Bakun Lake. The lake area variability is derived from 584 the high-resolution Landsat imagery data. The effective lake water level is determined using GRACE 585 data through a developed forward model. The land subsidence is then estimated by incorporating 586 Landsat and GRACE strengths into the developed disk load model, and the results are consistent with 587 the GPS measurement. 588 Recovering the local water storage is difficult using the GRACE data alone. This study demonstrates 589 how the low-resolution GRACE-derived ΔTWS data can be successfully used to determine and model 590 the surface water load over an area considerably less than the GRACE's intrinsic resolution, in 591 conjunction with high-resolution Landsat image data. Our approach provides the opportunity to 592 extract the local hydrological or geophysical signals that are still hidden in the GRACE observation,                  Spatial and temporal limited to a few hundred km and one month. Only provide the integrated water column, and cannot be disaggregated.

Landsat
Capture high spatial resolution surface reflectance, which can be effectively used to delineate the surface water area.
Only provide surface water component. Limited number of "cloud free" scenes.