Simulation of the Dynamic Water Storage and Its Gravitational E ﬀ ect in the Head Region of Three Gorges Reservoir Using Imageries of Gaofen-1

: The construction of a high-resolution dynamic water storage model, driven by the mass load of the huge water storage of the Three Gorges Reservoir (TGR), is the necessary basic data for accurately simulating changes in the geophysical ﬁeld, e.g., gravity, crustal deformation, and stress. However, previously established models cannot meet the needs of accurately simulating the impoundment e ﬀ ects of TGR, because these models were simpliﬁed and approximated and did not consider the variation of river boundaries caused by water level changes. In this study, we combined high-resolution Gaofen-1 (GF-1) satellite imageries and real-time water level in front of the dam and extracted 31 river boundaries of the head region of TGR between the lowest (145 m) and the highest (175 m) impoundment stages based on the Normalized Di ﬀ erential Water Index (NDWI) and threshold segmentation from Otsu method. Developed dynamic water storage model based on higher-resolution GF-1 data can show the true river boundary changes more exactly, especially in local areas. Compared to the previous approximate models, the model that we constructed accurately depicts the boundary distribution information of the di ﬀ erent impoundment stages. Moreover, we simulated TGR-induced gravitational e ﬀ ects based on the high-precision forward modeling of the dynamic water storage model (i.e., considering changes of dynamic water area and water level). The theoretical modelled results are consistent with in situ gravity measurements with the di ﬀ erence mainly within 10 µ Gal. Our results indicate that water storage variations of TGR mainly a ﬀ ect the gravity ﬁeld response within 1000 m of the reservoir bank with its maximum amplitude up to several hundred µ Gal. The dynamic water storage and its simulation results of gravitational e ﬀ ects can e ﬀ ectively eliminate the impact of surface water load driven by the TGR under human control and greatly improve the signal-to-noise ratio of regional gravity observational data. Thus, this work will be beneﬁcial in the application of geophysical and geodetic monitoring aimed to comprehensively track the local and regional geological structural stability, e.g., artiﬁcial reservoir induced earthquake and landslide.


Introduction
Determining the Earth's surface hydrological processes of terrestrial water storage or humaninduced water storage is considered a recent interest because the processes reflect the information of modern climate cycles changes [1]. The Three Gorges Reservoir (TGR) in China, covering the section from Yichang to Chongqing of Yangtze River basin (Figure 1), is one of the largest reservoirs in the world and a typical representative of man-made terrestrial water circulation [2,3]. Since the impoundment successfully reached its largest capacity (i.e., ~175 m water level) in 2008, the water level of the reservoir is perennially controlled in the periodic fluctuation between 145 m and 175 m for meeting the needs of flood control and irrigation in its lower reaches of the Yangtze River basin. The huge water storage change caused by impoundment in the TGR and its impacts on the regional environment [4][5][6] and accompanying variations in crustal load and deformation due to the water level rising and falling have been the focus issues for a long time [7,8]. The result of water storage fluctuation is considered as important inducing factors for geological hazards such as The huge water storage change caused by impoundment in the TGR and its impacts on the regional environment [4][5][6] and accompanying variations in crustal load and deformation due to the water level rising and falling have been the focus issues for a long time [7,8]. The result of water storage fluctuation is considered as important inducing factors for geological hazards such as reservoir earthquakes and landslides around the TGR [7,8]. Therefore, constructing a water storage model of the TGR is necessary to estimate the impacts due to water storage change. Models established previously [9][10][11][12] were relatively simple, and the variation of river boundaries caused by the impoundment was not considered. Since the change of reservoir water level has great influence on the nearshore area, where series of preprocessing before application. For removing the geometrical distortions, we applied the ASTER GDEM to orthorectify selected imageries. ASTER GDEM is a newer global DEM and has a ground resolution of ~30 m (https://lpdaac.usgs.gov/products/astgtmv002/). The imagery data preprocessing procedure shown in Figure 2 was completed using the ArcMap software (available at https://desktop.arcgis.com/zh-cn/arcmap/).   Satellite imagery compatible to our task needs to be selected from the large dataset of GF-1. The water level of the TGR varies throughout the whole year, but it usually stays at high water level for power generation in the winter months (November to February) and stays at low water level for flood control (June to August) [4]. Imageries acquired in these two periods are acceptable, except that there are a lot of clouds over the study area. Imageries acquired in the periods in which the water level quickly rises (September to October) and quickly falls (March to May) are not acceptable because those imageries cannot provide clear water boundaries due to the water level varying rapidly. Therefore, time frame for imagery acquisition is crucial in the study. We finally selected 11 imageries over the study area with little cloud cover in the period when water levels were about 145 m and 175 m ( Table 2).

Extraction of Water Boundaries
Water reflectance is generally rather low (4-5%) in the visible light range. Even though it is relatively high in the blue-green band, water reflectance is lower in the near infrared band and has a downtrend with increase in light wavelength. Hence, water bodies are mostly dark in the remote sensing imageries, and it is useful to distinguish the water body from other ground objects in the spectrum. In order to classify water targets from remote sense imageries, McFeeters [30] proposed the Normalized Differential Water Index (NDWI) applied to identify water bodies according to the difference of normalized ratio index of the green band and the near-infrared band. This approach has been previously used for distinguishing and extracting the water boundaries [22,[31][32][33]. We use NDWI [30] to distinguish water bodies by applying the following equation: where B2 and B4 denote reflectance in blue-green band and near-infrared band, respectively. NDWI can enhance the difference between B2 and B4 and effectively separate water from other land covers, and thereby, it can be taken as an indicator for discriminating water and non-water features. The water boundary extraction of the corrected fused multispectral imagery can be implemented by the following steps ( Figure 3): (1) calculating the NDWI imagery, (2) using the Otsu method [34] to get the optimal threshold and segment the water body information, (3) through the local visual interpretation adjustment, we can get the corrected water boundaries.
Remote Sens. 2020, x, x FOR PEER REVIEW 7 of 23  Figure 4 shows a part of the fused multispectral imagery (Figure 4a) and the NDWI imagery ( Figure 4b). It is further necessary for selecting an appropriate threshold to distinguish water from other objects in NDWI imagery. The Otsu method [34] is considered as an efficient approach to perform binarized segmentation that is usually applied to identify water bodies or non-water bodies  Figure 4b). It is further necessary for selecting an appropriate threshold to distinguish water from other objects in NDWI imagery. The Otsu method [34] is considered as an efficient approach to perform binarized segmentation that is usually applied to identify water bodies or non-water bodies [22].
Remote Sens. 2020, x, x FOR PEER REVIEW 6 of 20 perform binarized segmentation that is usually applied to identify water bodies or non-water bodies [22]. The threshold value is calculated based on the rule of maximum between-class variance and minimum within-class variance in the Otsu method. Following this principle, assume [a, b] as the range of index value in NDWI imagery; any threshold t (a ≤ t ≤ b) may divide indices into two categories, i.e., indices are attributed to the water category as they are larger than t, and others are attributed to the non-water category. The optimal threshold should be calculated in equations by following: where σt is the between-class variance with respect to t; Pw (water) and Pnw (non-water) are probabilities of pixel units divided into two categories respectively, and Mw (water) and Mnw (nonwater) are their average NDWI values. Given different t for test, the optimal threshold Topt can be determined by: The spectrum divergence among different imageries is usually inconsistent due to different observation periods and locations. Thus, selecting an optimal threshold for each imagery is required. We change threshold t with an interval of 0.01 and iteratively search maximum between-class variances for every imagery. Achieved optimal thresholds for imageries used in this study are listed in Table 3. The threshold value is calculated based on the rule of maximum between-class variance and minimum within-class variance in the Otsu method. Following this principle, assume [a, b] as the range of index value in NDWI imagery; any threshold t (a ≤ t ≤ b) may divide indices into two categories, i.e., indices are attributed to the water category as they are larger than t, and others are attributed to the non-water category. The optimal threshold should be calculated in equations by following: where σ t is the between-class variance with respect to t; P w (water) and P nw (non-water) are probabilities of pixel units divided into two categories respectively, and M w (water) and M nw (non-water) are their average NDWI values. Given different t for test, the optimal threshold T opt can be determined by: The spectrum divergence among different imageries is usually inconsistent due to different observation periods and locations. Thus, selecting an optimal threshold for each imagery is required. We change threshold t with an interval of 0.01 and iteratively search maximum between-class variances for every imagery. Achieved optimal thresholds for imageries used in this study are listed in Table 3. Binarized imagery, according to the achieved optimal threshold, can be obtained by applying the following formula: where A w denotes the attribution value of water in the binarized segmentation imagery. Water boundaries at pixel level can be easily extracted from a binarized imagery through the conversion from raster imagery to vector map. Boundary lines are composed of segments that connect nodes between adjacent water and non-water pixels, in which vector boundary line dataset consists of node coordinates. Figure 4c shows a part of binarized imagery transformed from a NDWI imagery in this study.
Although algorithms for the automatic target identification are efficient and successful, they are often influenced by interferences such as mountain shadows, ships, bridges, and floating objects ( Figure 4c). Accordingly, further corrections should be manually accomplished based on the visual inspection of the imageries as a complement [35][36][37] (Figure 4d). The extracted water boundary lines at water levels of about 145 m and 175 m are presented in Figure 5.

Dynamic Model Construction
Previously applied methods [9][10][11][12] for estimating the difference of water storage in the TGR measured the variation of water volume in which the product between the water surface area of closed water boundary lines (at the lowest or the highest), and the water level difference was taken as the estimation of volume variation. However, riverbank and riverbed topography in the TGR area was not considered. The original terrain of the TGR area, before impoundment, was rough and the details can be seen in Figure 6a. It can be noticed that several underwater hills are not fully submerged during the highest and lowest water levels, while others are only submerged during the highest water levels.
Remote Sens. 2020, x, x FOR PEER REVIEW 8 of 20

Dynamic Model Construction
Previously applied methods [9][10][11][12] for estimating the difference of water storage in the TGR measured the variation of water volume in which the product between the water surface area of closed water boundary lines (at the lowest or the highest), and the water level difference was taken as the estimation of volume variation. However, riverbank and riverbed topography in the TGR area was not considered. The original terrain of the TGR area, before impoundment, was rough and the details can be seen in Figure 6a. It can be noticed that several underwater hills are not fully submerged during the highest and lowest water levels, while others are only submerged during the highest water levels. As hills are only flooded at the highest water level (Figure 6a, situation II), the water storage may be miscounted. Acquiring other water boundary lines between the highest and the lowest water level can be helpful in diminishing these impacts and miscounts. Since satellite imageries are unsuitable to acquire water boundary lines during water level rising and falling in the TGR, the interpolation based on the lowest and highest water levels may be a reasonable approach. We simulated the bank terrain between water boundary lines at water levels of 145 m and 175 m by the minimum curvature algorithm in order to consider the topography of the reservoir banks. Intermediate water boundaries were obtained from the simulated surface of banks. Afterwards, the water boundary lines with a uniform interval between the highest and the lowest water levels could be acquired. Figure 6b shows the example near the front of the Three Gorges Dam. The results show that the highest and the lowest water level boundaries well constrain the distribution of the intermediate water level lines. Finally, a reasonable water storage model is obtained and it reflects the dynamic process of gradual changes in water level.
Meticulously divided water boundaries provide more possibilities to construct a detailed water storage model at different water levels in the TGR. According to the obtained water boundary lines, the water body of the TGR was divided into a series of horizontal layers. The partial or total combination of layers represents the water storage at different water levels, and it can be used for the estimation of water load and other geophysical effects at various water levels.
Geometrically, each layer of the water body can be described as a complex shape body closed by two respective horizontal planes on the top and bottom and rough surfaces on two reservoir banks, or it can be defined as a polyhedron. For comprehending the construction of a layer of the water body, we detached the river into a series of small segments.
The layered water body shaped in a polyhedron (e.g., the single-layer water body in Figure 7a represents one layer of the multi-layer model such as Layer 3 in Figure 7b) can be separated into three parts: one quadrilateral prism and two pentahedrons (Figure 7c). The quadrilateral prism is the major As hills are only flooded at the highest water level (Figure 6a, situation II), the water storage may be miscounted. Acquiring other water boundary lines between the highest and the lowest water level can be helpful in diminishing these impacts and miscounts. Since satellite imageries are unsuitable to acquire water boundary lines during water level rising and falling in the TGR, the interpolation based on the lowest and highest water levels may be a reasonable approach. We simulated the bank terrain between water boundary lines at water levels of 145 m and 175 m by the minimum curvature algorithm in order to consider the topography of the reservoir banks. Intermediate water boundaries were obtained from the simulated surface of banks. Afterwards, the water boundary lines with a uniform interval between the highest and the lowest water levels could be acquired. Figure 6b shows the example near the front of the Three Gorges Dam. The results show that the highest and the lowest water level boundaries well constrain the distribution of the intermediate water level lines. Finally, a reasonable water storage model is obtained and it reflects the dynamic process of gradual changes in water level.
Meticulously divided water boundaries provide more possibilities to construct a detailed water storage model at different water levels in the TGR. According to the obtained water boundary lines, the water body of the TGR was divided into a series of horizontal layers. The partial or total combination of layers represents the water storage at different water levels, and it can be used for the estimation of water load and other geophysical effects at various water levels.
Geometrically, each layer of the water body can be described as a complex shape body closed by two respective horizontal planes on the top and bottom and rough surfaces on two reservoir banks, or it can be defined as a polyhedron. For comprehending the construction of a layer of the water body, we detached the river into a series of small segments.
The layered water body shaped in a polyhedron (e.g., the single-layer water body in Figure 7a represents one layer of the multi-layer model such as Layer 3 in Figure 7b) can be separated into three parts: one quadrilateral prism and two pentahedrons (Figure 7c). The quadrilateral prism is the major part of the layer in the segment. After connecting neighboring quadrilateral prisms in segments, those quadrilateral prisms are combined to a huge polygonal prism, with the consistent area on the top and bottom and a certain thickness along the total river (see the right side in Figure 7d). This is a relatively simple model that allows convenient calculation of its load and gravitational effect. However, this model is not appropriate for total water storage calculation. In that case, pentahedrons on the two reservoir banks need to be considered. For example, each considered pentahedron shape is represented by closed surface with one horizontal plane on the top, three vertical planes on the sides, and one quadrangle facet on the bank.
Remote Sens. 2020, x, x FOR PEER REVIEW 9 of 20 part of the layer in the segment. After connecting neighboring quadrilateral prisms in segments, those quadrilateral prisms are combined to a huge polygonal prism, with the consistent area on the top and bottom and a certain thickness along the total river (see the right side in Figure 7d). This is a relatively simple model that allows convenient calculation of its load and gravitational effect. However, this model is not appropriate for total water storage calculation. In that case, pentahedrons on the two reservoir banks need to be considered. For example, each considered pentahedron shape is represented by closed surface with one horizontal plane on the top, three vertical planes on the sides, and one quadrangle facet on the bank. The volume depends on the bank terrain and the corners (vertexes) of the facet (AHG'B' or DEF'C' in Figure 7c) are perhaps non-coplanar, i.e., the body is not an exact pentahedron. By subdividing the facet AHG'B' into two triangular facets (ΔAB'H and ΔHB'G'), the pentahedron becomes an asymmetric hexahedron or polyhedron, and this assumption is suitable for all these bodies on the bank in each segment. By connecting the neighbor asymmetric hexahedrons in all segments one by one, they will assemble two hexahedron or polyhedron chains along reservoir banks in which the immerged bank surface may be subdivided into planar triangle facets. Bank surface between the adjacent two water boundary lines can be simulated by applying the Delaunay triangulation along reservoir banks, e.g., Si is the ith triangular facet, and Pi,1, Pi,2 and Pi,3 are the corners of the facet defined by the Delaunay triangulation algorithm (see left part of Figure 7d). Finally, the water storage at arbitrary water level in the TGR is the summary of volume from the lowest layer to the layer at the correspondent water level. This approach is helpful for accurately investigating the variation of the water storage and its geophysical responses in the TGR. The volume depends on the bank terrain and the corners (vertexes) of the facet (AHG'B' or DEF'C' in Figure 7c) are perhaps non-coplanar, i.e., the body is not an exact pentahedron. By subdividing the facet AHG'B' into two triangular facets (∆AB'H and ∆HB'G'), the pentahedron becomes an asymmetric hexahedron or polyhedron, and this assumption is suitable for all these bodies on the bank in each segment. By connecting the neighbor asymmetric hexahedrons in all segments one by one, they will assemble two hexahedron or polyhedron chains along reservoir banks in which the immerged bank surface may be subdivided into planar triangle facets. Bank surface between the adjacent two water boundary lines can be simulated by applying the Delaunay triangulation along reservoir banks, e.g., S i is the ith triangular facet, and P i,1 , P i,2 and P i,3 are the corners of the facet defined by the Delaunay triangulation algorithm (see left part of Figure 7d). Finally, the water storage at arbitrary water level in the TGR is the summary of volume from the lowest layer to the layer at the correspondent water level. This approach is helpful for accurately investigating the variation of the water storage and its geophysical responses in the TGR.

Principle of Gravitational Effect Computation
As discussed before, the layered polyhedron is the basic element for constructing the dynamic model of the water storage in the TGR. Previous studies [38,39] provided completely analytical expressions for the first and second derivatives of the gravitational potential in arbitrary directions generated by a homogeneous polyhedral body closed by polygonal facets.
Usually, the gravitational potential U, generated by a homogeneous body, can be given in integral form as where G is the constant of universal gravitation, ρ is the density of the body, V is the volume of the body, and r is the distance from origin point to integral element central point (x, y, z), i.e., r = (x 2 + y 2 + z 2 ) 1/2 . The derivative of gravitational potential in arbitrary directions is the component of the attraction force in that direction. If U k denotes the component of the attraction force in k-direction, it can be expressed as where the sign ∇ denotes gradient operator, k is the unit vector in k-direction. Applying the Gauss's divergence theorem, Equation (6) becomes where S is the closed surface of body V, n is the unit vector in outward normal direction on S. The S composes of finite number polygonal facets. Equation (7) can be rewritten in: where S i is the ith facet, and n i is the unit vector in outward normal direction on S i . For solving the surface integral in Equation (8), the coordinates may rotate to another Cartesian system (X, Y, Z), in which the Z-direction is coincident with n i . The surface integral in the XOY plane is transformed to line integral along the closed edge of S i based on the Green Formula. Because the edge composes of finite line-segments, line integral can be carried out by subsection integral. Hence, the surface integral of abbreviated expression is where J j (i) indicates the line integral from jth vertex to (j + 1)th vertex on the edge of S i . Vertexes must be numbered clockwise [39], and for the last edge, vertex (j + 1) corresponds to vertex 1. See Appendix A for detailed deducing process. Finally, Equation (8) can be expressed as the summary of line integrals by following: Remote Sens. 2020, 12, 3353 11 of 20 The gravitational effect (g) usually represents the downward vertical component of the attraction force and the Equation (10) can be rewritten as where φ i is the spatial angle between the outward normal direction of S i and downward vertical direction. The Equation (11) is appropriate for computing the gravitational effects of the whole layered water body that includes the middle part of polygonal prism and polyhedron chains on the both banks, or only a single polyhedron. The gravitational effect of the water body at a certain water level can be obtained by accumulating those from the lowest water level.
Two simulating approaches are applied in this study. A synthetic model is designed and shown in Figure 8 in order to manifest the difference of the gravitational effects generated by water bodies. An asymmetry heptahedron (composed of a cuboid and an asymmetry hexahedron on the left, with regard to rough bank surface) is applied to simulate the model. The model density is defined as 1.0 × 10 3 kg/m 3 (i.e., fresh water). We selected observational points P 1 (0,0,0), P 2 (0,500,0), and P 3 (0,1000,0), and the distance unit is meters. Afterwards, the gravitational effects at these points are calculated and listed in Table 4. These results show that using the cuboid for simulation and without considering water on reservoir bank would result in loss of approximately 47% in gravitational effect near the riverside. Therefore, for comprehensive estimation of the gravitational effect, it is necessary to use polyhedrons to simulate water body instead of prism bodies.
Remote Sens. 2020, x, x FOR PEER REVIEW 11 of 20 The Equation (11) is appropriate for computing the gravitational effects of the whole layered water body that includes the middle part of polygonal prism and polyhedron chains on the both banks, or only a single polyhedron. The gravitational effect of the water body at a certain water level can be obtained by accumulating those from the lowest water level.
Two simulating approaches are applied in this study. A synthetic model is designed and shown in Figure 8 in order to manifest the difference of the gravitational effects generated by water bodies. An asymmetry heptahedron (composed of a cuboid and an asymmetry hexahedron on the left, with regard to rough bank surface) is applied to simulate the model. The model density is defined as 1.0 × 10 3 kg/m 3 (i.e., fresh water). We selected observational points P1(0,0,0), P2(0,500,0), and P3(0,1000,0), and the distance unit is meters. Afterwards, the gravitational effects at these points are calculated and listed in Table 4. These results show that using the cuboid for simulation and without considering water on reservoir bank would result in loss of approximately 47% in gravitational effect near the riverside. Therefore, for comprehensive estimation of the gravitational effect, it is necessary to use polyhedrons to simulate water body instead of prism bodies.

Gravitational Effects of Dynamic Water Storage in Head Region of TGR
Using the approach proposed in Section 2.4, we calculated the gravitational effects at different water impoundment periods based on the obtained water boundaries. In the theoretic modeling, the lowest and the highest water levels are defined at 145 m and 175 m, respectively. The gravitational effects on the surface in the head region of TGR are computed, for which SRTM DEM is employed to approximate the surface with a ground resolution of about 30 m. Figure 9 shows results during six water impoundment periods. The gravitational effects increase with rising water levels, and the maximum of more than 600 µGal is obtained at period 145-175 m. Gravitational effects are negative at a few places because the elevation of observation points is lower than the water level at that place. There is a tendency for high gravitational effects to appear near the reservoir, and their intensity rapidly decays as distance from the reservoir bank increases. Furthermore, the impact of local water storage on gravitational effect is significantly increased in areas where tributaries meet. In general, the gravitational effects of water storage are affected by the relative distance of the local reservoir bank and the relative elevation difference between the calculation site and water level surface.
Remote Sens. 2020, x, x FOR PEER REVIEW 12 of 20 effects on the surface in the head region of TGR are computed, for which SRTM DEM is employed to approximate the surface with a ground resolution of about 30 m. Figure 9 shows results during six water impoundment periods. The gravitational effects increase with rising water levels, and the maximum of more than 600 μGal is obtained at period 145-175 m. Gravitational effects are negative at a few places because the elevation of observation points is lower than the water level at that place. There is a tendency for high gravitational effects to appear near the reservoir, and their intensity rapidly decays as distance from the reservoir bank increases. Furthermore, the impact of local water storage on gravitational effect is significantly increased in areas where tributaries meet. In general, the gravitational effects of water storage are affected by the relative distance of the local reservoir bank and the relative elevation difference between the calculation site and water level surface.

Comparison with the Measurement of Absolute Gravimeter
We selected four field gravity observation sites in the head region of TGR ( Figure 5) that can be used for detecting regional gravity changes caused by impoundment. Afterwards, we calculated the time-varying gravitational effects based on dynamic water storage model at in situ observation sites ( Figure 10). Since 2010, field observations of absolute gravity were carried out twice per year, i.e., in summer (low water level, ~145 m) and in winter (high water level, ~175 m). The field observations included vertical gravity gradient measurements from a drop cart position to the ground. This study used repeated measurements of absolute gravity to verify the validity and accuracy of forward

Comparison with the Measurement of Absolute Gravimeter
We selected four field gravity observation sites in the head region of TGR ( Figure 5) that can be used for detecting regional gravity changes caused by impoundment. Afterwards, we calculated the time-varying gravitational effects based on dynamic water storage model at in situ observation sites ( Figure 10). Since 2010, field observations of absolute gravity were carried out twice per year, i.e., in summer (low water level,~145 m) and in winter (high water level,~175 m). The field observations included vertical gravity gradient measurements from a drop cart position to the ground. This study used repeated measurements of absolute gravity to verify the validity and accuracy of forward modeling from the dynamic water storage model. The results show that the A10 observations display a similarly regular cyclic gravitational effect due to periodical water storage changes controlled by the dam (Figure 10). Furthermore, the theoretical gravitational effects calculated by the dynamic water storage model are highly consistent with in situ measured A10 data. It should be noted that the measurements of A10 represent the absolute gravity values, while the calculated gravitational effect of the water body above the lowest water level is a relative value obtained by the dynamic water storage model. Thus, we have separately calculated the average values of the measured absolute gravity during the lowest and the highest water level impoundment period. For example, for the A10 observation data acquired during the lowest water level stage, we used the observed absolute gravity data (gobserved in 145m) minus the average value of all observed data (gobserved mean in 145m) during the lowest water level stage; then, we added the average value of the theoretical gravitational effects (gcalculated mean in 145m) in the corresponding observation time, as a comparison reference value of the measured data (i.e., gobserved_lowest = gobserved in 145m − gobserved mean in 145m + gcalculated mean in 145m). Afterwards, we used a similar process for the observation data acquired during the highest water level stage (i.e., gobserved_highest = gobserved in 175m − gobserved mean in 175m + gcalculated mean in 175m). The results show that the amplitude of gravitational effects varies greatly at different sites, which mainly depends on the mass changes and the relative position between A10 sites and the dynamic water storage model. The observed absolute gravity data verified the simulation predicted results with small deviations of a few to dozens of μGal. It has to be noted that the A10 observed results also include other potential environmental effects (e.g., climate-induced rainfall, groundwater storage variations, etc.), in addition to the direct gravitational effects of TGR impoundment and A10 instrument measurement uncertainty (10 μGal). These influencing factors will be discussed in the next section. The results show that the A10 observations display a similarly regular cyclic gravitational effect due to periodical water storage changes controlled by the dam (Figure 10). Furthermore, the theoretical gravitational effects calculated by the dynamic water storage model are highly consistent with in situ measured A10 data. It should be noted that the measurements of A10 represent the absolute gravity values, while the calculated gravitational effect of the water body above the lowest water level is a relative value obtained by the dynamic water storage model. Thus, we have separately calculated the average values of the measured absolute gravity during the lowest and the highest water level impoundment period. For example, for the A10 observation data acquired during the lowest water level stage, we used the observed absolute gravity data (g observed in 145m ) minus the average value of all observed data (g observed mean in 145m ) during the lowest water level stage; then, we added the average value of the theoretical gravitational effects (g calculated mean in 145m ) in the corresponding observation time, as a comparison reference value of the measured data (i.e., g observed_lowest = g observed in 145m − g observed mean in 145m + g calculated mean in 145m ). Afterwards, we used a similar process for the observation data acquired during the highest water level stage (i.e., g observed_highest = g observed in 175m − g observed mean in 175m + g calculated mean in 175m ). The results show that the amplitude of gravitational effects varies greatly at different sites, which mainly depends on the mass changes and the relative position between A10 sites and the dynamic water storage model. The observed absolute gravity data verified the simulation predicted results with small deviations of a few to dozens of µGal. It has to be noted that the A10 observed results also include other potential environmental effects (e.g., climate-induced rainfall, groundwater storage variations, etc.), in addition to the direct gravitational effects of TGR impoundment and A10 instrument measurement uncertainty (10 µGal). These influencing factors will be discussed in the next section.

The Capacity of the Dynamic Water Storage Model
The dynamic water storage model constructed in the head region of TGR is mainly composed of the 31 water boundaries of the reservoir and the real-time water level in front of the dam. Accordingly, this model is based on the valid distribution pattern of the water body changes due to the artificially controlled dam. Compared with the previous static model constructed by SRTM DEM data (about 90 m resolution), our dynamic model is based on the GF-1 imagery that represents the changes of the water boundaries (more fit with the real water body) in different periods and with better resolution. For the variation of water storage in the head region of TGR, we calculated the volumes and inundated areas at different water levels by applying our dynamic water storage model and previous static model based on SRTM DEM (Figure 11).  The calculated volumes from two different water storage models are similar, yet the variation characteristics of inundated areas are different. The horizontal changes in riverbank lines basically vary on a few dozens of meters ( Figure 5). Accordingly, the lower resolution of SRTM DEM (~90 m) in the previous static model can produce errors when recognizing river boundaries. This approximate result based on STRM may not show significant differences when we use the entire TGR hydrological information; however, it will show obvious impact in a small local area such as the head region of TGR. For example, if we assume that the water storage model adds 1 grid after the water level is increased by 1 m, then the increased area is mainly based on the horizontal resolution of the grid, i.e., the grid with 1-m resolution will cause an increase in the area of 1 m × 1 m, while a grid with a STRM 90-m resolution will result in an increase of 90 m × 90 m. For the local area selected in this study, our dynamic water storage model based on higher-resolution GF-1 data can provide more detailed river boundary changes.

Simulation Accuracy
The gravitational effect of a water storage model depends on the characteristics of the model itself and the spatial position between the observation site and the model, i.e., the difference in The calculated volumes from two different water storage models are similar, yet the variation characteristics of inundated areas are different. The horizontal changes in riverbank lines basically vary on a few dozens of meters ( Figure 5). Accordingly, the lower resolution of SRTM DEM (~90 m) in the previous static model can produce errors when recognizing river boundaries. This approximate result based on STRM may not show significant differences when we use the entire TGR hydrological information; however, it will show obvious impact in a small local area such as the head region of TGR. For example, if we assume that the water storage model adds 1 grid after the water level is increased by 1 m, then the increased area is mainly based on the horizontal resolution of the grid, i.e., the grid with 1-m resolution will cause an increase in the area of 1 m × 1 m, while a grid with a STRM 90-m resolution will result in an increase of 90 m × 90 m. For the local area selected in this study, our dynamic water storage model based on higher-resolution GF-1 data can provide more detailed river boundary changes.

Simulation Accuracy
The gravitational effect of a water storage model depends on the characteristics of the model itself and the spatial position between the observation site and the model, i.e., the difference in distance and height between the calculation point and the model. In this study, we tested the impacts of the forward modeling accuracy due to different resolution of model boundary. Firstly, we applied a dynamic polygonal-slab model that only uses the lowest and the highest water boundaries as the model boundary. This simplified slab model can approximate slope surface to vertical plane when calculating the gravitational effect, yet it neglects the influence of the millions of triangular facets of the slope (see Figure 7d). Afterwards, we calculated and compared the gravitational effects in the surrounding region based on the dynamic polygonal-slab model (Figure 12a) and the precise dynamic water storage model (Figure 12b). The comparison results show that the dynamic polygonal-slab model will cause a large calculation error when the water boundary near the reservoir bank drastically changes (Figure 12c). In other words, when the area is closer to the water body, there is also larger amplitude in gravitational effect changes and the calculation error in simplified models. This is especially evident in the tributary area where the calculation error changes faster. The dynamic polygonal-slab model is characterized by the fast calculation speed and is suitable for the analysis of the distant area from the river bank (e.g., distance over 2000 m), where the amplitude of gravitational effect of the water storage model is relatively small and slowly fluctuates (Figure 9h).
Remote Sens. 2020, x, x FOR PEER REVIEW 15 of 20 surrounding region based on the dynamic polygonal-slab model (Figure 12a) and the precise dynamic water storage model (Figure 12b). The comparison results show that the dynamic polygonal-slab model will cause a large calculation error when the water boundary near the reservoir bank drastically changes (Figure 12c). In other words, when the area is closer to the water body, there is also larger amplitude in gravitational effect changes and the calculation error in simplified models. This is especially evident in the tributary area where the calculation error changes faster. The dynamic polygonal-slab model is characterized by the fast calculation speed and is suitable for the analysis of the distant area from the river bank (e.g., distance over 2000 m), where the amplitude of gravitational effect of the water storage model is relatively small and slowly fluctuates (Figure 9h). For the situation at the single observation site (e.g., site 5 in Figure 5), we calculated the timevarying gravitational effects based on the precise dynamic water storage model and the dynamic polygonal-slab model, respectively (Figure 12d). The distance of horizontal position and elevation difference between site 5 and the nearest model boundary ( Figure 5) are ~349 m and ~93 m (relative to 145 m), respectively. The maximum variation amplitude of the theoretical gravitational effect at site 5 can reach up to 70 μGal based on the dynamic water storage model, and the maximum difference can reach up to 4.58 μGal. It is noticed that the difference between the two water storage models will be larger when the observation site or the calculation point is closer to the river bank. Moreover, gravitational effects based on two water model appear to be significantly inconsistent with different water levels, e.g., the difference between the forward modeling results of the two models is relatively large during sharp fluctuations of water level (releasing water in spring and flooding in For the situation at the single observation site (e.g., site 5 in Figure 5), we calculated the time-varying gravitational effects based on the precise dynamic water storage model and the dynamic polygonal-slab model, respectively (Figure 12d). The distance of horizontal position and elevation difference between site 5 and the nearest model boundary ( Figure 5) are~349 m and~93 m (relative to 145 m), respectively. The maximum variation amplitude of the theoretical gravitational effect at site 5 can reach up to 70 µGal based on the dynamic water storage model, and the maximum difference can reach up to 4.58 µGal. It is noticed that the difference between the two water storage models will be larger when the observation site or the calculation point is closer to the river bank. Moreover, gravitational effects based on two water model appear to be significantly inconsistent with different water levels, e.g., the difference between the forward modeling results of the two models is relatively large during sharp fluctuations of water level (releasing water in spring and flooding in monsoon) and relatively small during stable water level periods. In this case, the precise dynamic water storage model can identify with greater accuracy the local boundary changes and is more suitable for precise correction of observational points around the reservoir banks in order to extract other hydrological or geophysical signals, e.g., groundwater leakage and tectonic movement.

Other Influences
We noticed that our model simulation results are in good agreement with the gravity observational data at larger gravity difference scales (e.g., dozens or even hundreds of µGal as shown in Figure 10). This indicates that the change in gravitational effects caused by the impoundment is the main factor for the non-tidal gravity changes within 1000 m of the reservoir bank, with a few subtle differences between modelled and observed results. In fact, A10 absolute gravity observations reflect all gravity effects including TGR water impoundment and other surface processes around the site, after performed preprocessing (e.g., removal of effects such as solid earth tides, ocean tides, pole motion, normal field, and elevation correction at each site). For example, several unknown environmental gravity changes such as groundwater changes [5], seasonal changes in rainfall [4], and tectonic movements [13] are reflected in A10 results, and they can cause gravity changes of about 10 µGal [28]. Wang et al. [5] showed that the interannual change of the equivalent groundwater height in the TGR area is about 1.6 m, which leads to the observed gravitational changes of up to tens of µGal.
The elastic deformation of the Earth caused by the huge impoundment water load can also cause the change of the gravitational field [12], yet its amplitude is relatively small (within 6 µGal). In addition, regional economic development (e.g., construction of roads, bridges, buildings, etc.) and local topographic changes caused by geological disasters like landslides are also factors that can cause disagreement between simulated and measured results. Overall, these factors should be considered in future study, e.g., by separating secondary signals directly or indirectly from the water storage changes in TGR.

Conclusions
In this study, we used GF-1 high-resolution satellite imageries and water level in front of the dam in order to develop a dynamic water storage model which would provide support for the fine simulation of the surface dynamics of TGR. By utilizing the forward modeling method of complex body, it was possible to obtain the high-precision simulation results of gravitational effects that provide important real-time hydrologic correction data for identifying and distinguishing the gravitational effects between TGR-induced water changes and other hydrological elements. The main conclusions are as follows: 1.
In front of the dam (Zigui area), the annual water boundary change can be up to 100 m in horizontal direction. In the tributary area of the TGR, the changes of the water boundary are more obvious, especially in the north of Badong area. The dynamic water storage model constructed in this study consists of water boundaries (31 in total) at different impoundment periods, i.e., from 145 m to 175 m. The resolution of the model is approximately 8 m in horizontal direction and approximately 1 m in vertical direction, which is an important way to describe the local dynamic characteristics in the impounding and releasing water process.

2.
The high-precision forward modeling method of gravitational effect can avoid the error caused by model approximation. The gravitational effect of TGR water storage is obvious within 1000 m of the reservoir bank, and the maximum gravity change caused by the 30-m water level change (from 145 m to 175 m) can reach nearly 600 µGal. The ground multiyear gravity observation data from A10 absolute gravimeters verified the simulations. Our results confirm that the gravitational effects change caused by TGR impoundment is the main contribution of the gravity field changes in the study area. The differences between in situ observations and simulation results can help in further tracking of the effects of geological disasters, such as earthquakes and landslides. 3.
The comparison results show that the dynamic water storage model has more advantages in identifying changes in water storage in local areas and is statistically more accurate in obtaining TGR capacity information when compared to the previous model. However, the potential users should consider the spatial position and coordinate accuracy of the calculated point when applying the dynamic water storage model. 4.
The accuracy of the dynamic water storage model is related to factors such as satellite imagery preprocessing (e.g., data calibration), imagery quality, cloud cover, and extraction method for water boundaries. In addition, manual visual interpretation can solve the influence of interference factors (e.g., river vessels, mountain shadows, water pollutants, etc.) With the continuous accumulation and update of the remote sensing data and the development of water boundary extraction technology, it can be expected that the accuracy of the water storage model will be improved in the future.  Before discussing the definite integral in ζ-axis, we have to define the function V This V can easily be rewritten as Thus, the line integral J j (i) in Equation (A2) can be expressed as J j (i) = − cos ψ ξ j+1 ξ j Vdξ = (ξ cos ψ − η sin ψ) ln[ξ sin ψ + η cos ψ + (ξ 2 + η 2 + Z 2 ) 1/2 ] − ξ cos ψ + η ln[ξ + (ξ 2 + η 2 + Z 2 ) 1/2 ] Finally, the line integral J j (i) of abbreviated expression is J j (i) = (−X sin ψ + Y sin ψ) ln[X cos ψ + Y sin ψ + (X 2 + Y 2 + Z 2 ) 1/2 ] +2Z tan −1 (1 + sin ψ)[Y + (X 2 + Y 2 + Z 2 ) 1/2 ] + X cos ψ Z cos ψ