remote Assessing Surface Water Losses and Gains under Rapid Urbanization for SDG 6.6.1 Using Long-Term Landsat Imagery in the Guangdong-Hong Kong-Macao Greater Bay Area, China

: As one of the most open and dynamic regions in China, the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) has been urbanizing rapidly in recent decades. The surface water in the GBA also has been suffering from urbanization and intensiﬁed human activities. The study aimed to characterize the spatiotemporal patterns and assess the losses and gains of surface water caused by urbanization in the GBA via long time-series remote sensing data, which could support the progress towards sustainable development goals (SDGs) set by the United Nations, especially for measuring SDG 6.6.1 indicator. Firstly, utilizing 4750 continuous Landsat TM/ETM+/OLI images during 1986–2020 and the Google Earth Engine cloud platform, the multiple index water detection rule (MIWDR) was performed to extract surface water extent in the GBA. Secondly, we achieved surface water dynamic type classiﬁcation based on annual water inundation frequency time-series in the GBA. Finally, the spatial distribution and temporal variation of urbanization-induced water losses and gains were analyzed through a land cover transfer matrix. Results showed that (1) the average minimal and maximal surface water extents of the GBA during 1986–2020 were 2017.62 km 2 and 6129.55 km 2 , respectively. The maximal surface water extent fell rapidly from 7897.96 km 2 in 2001 to 5087.46 km 2 in 2020, with a loss speed of 155.41 km 2 per year ( R 2 = 0.86). (2) The surface water areas of permanent and dynamic types were 1529.02 km 2 and 2064.99 km 2 during 2000–2020, accounting for 42.54% and 57.46% of all water-related areas, respectively. (3) The surface water extent occupied by impervious land surfaces showed a signiﬁcant linear downward trend ( R 2 = 0.98, slope = 36.41 km 2 per year), while the surface water restored from impervious land surfaces denoted a slight growing trend ( R 2 = 0.86, slope = 0.99 km 2 per year). Our study monitored the long-term changes in the surface water of the GBA, which can provide valuable information for the sustainable development of the GBA urban agglomeration. In addition, the proposed framework can easily be implemented in other similar regions worldwide.


Introduction
Terrestrial surface water bodies, which include lakes, rivers, reservoirs, ponds and partially impounded farmland, are a critically important component of water resources for human life and terrestrial ecosystems [1]. Surface water provides and supports diverse ecosystem services for human activities and natural habitats [2].
(2) to achieve surface water dynamic type classification and mapping based on multiyear water inundation frequency time-series; (3) to analyze the impact of urbanization on surface water losses and gains quantitively from the aspect of realizing SDG 6.6.1 and discuss the possible climate causes for water area changes in the GBA.

Study Area
Located in the southern coastal area of China, the GBA (latitude: 111 • 20 E-115 • 24 E, longitude: 21 • 32 N-24 • 26 N) is a region across nine cities of Guangdong Province and two special administrative regions, which includes Zhaoqing, Guangzhou, Huizhou, Foshan, Dongguan, Jiangmen, Zhongshan, Shenzhen, Hong Kong and Macao, as shown in Figure 1. The total area of the GBA is approximately 56,000 km 2 . In addition, the GBA is situated in the humid subtropical monsoon climate zone, with abundant rainfall during rainy and hot periods (from April to September); its annual precipitation is approximately 1300~2500 mm, and its annual average temperature is approximately 22 • C [34][35][36].
traction method to explore its spatial variations and temporal change trends (by OLS-RS method) of surface water extent in the GBA; (2) to achieve surface water dynamic type classification and mapping based on multiyear water inundation frequency time-series; (3) to analyze the impact of urbanization on surface water losses and gains quantitively from the aspect of realizing SDG 6.6.1 and discuss the possible climate causes for water area changes in the GBA.

Study Area
Located in the southern coastal area of China, the GBA (latitude: 111°20′ E-115°24′ E, longitude: 21°32′ N-24°26′ N) is a region across nine cities of Guangdong Province and two special administrative regions, which includes Zhaoqing, Guangzhou, Huizhou, Foshan, Dongguan, Jiangmen, Zhongshan, Shenzhen, Hong Kong and Macao, as shown in Figure 1. The total area of the GBA is approximately 56,000 km 2 . In addition, the GBA is situated in the humid subtropical monsoon climate zone, with abundant rainfall during rainy and hot periods (from April to September); its annual precipitation is approximately 1300~2500 mm, and its annual average temperature is approximately 22 °C [34][35][36].
With a total population of more than 86 million people and regional GDP of approximately USD 1668.8 billion at the end of 2020, which was equivalent to 12% of the GDP for all of mainland China, the GBA is the largest and richest economic region in South China. The GBA is also the starting point of the Maritime Silk Road and the gateway of China to the world [37].  With a total population of more than 86 million people and regional GDP of approximately USD 1668.8 billion at the end of 2020, which was equivalent to 12% of the GDP for all of mainland China, the GBA is the largest and richest economic region in South China. The GBA is also the starting point of the Maritime Silk Road and the gateway of China to the world [37].

Landsat Data
In this study, we selected Landsat satellite images to extract water areas in the GBA. With a 30-m spatial resolution and 16-day temporal revisit period, all available Landsat 5 TM, 7 ETM+ and 8 OLI surface reflectance images covering the entire GBA from 1986 to 2020 were utilized and processed on the GEE platform for convenience and efficiency. The temporal distribution of image quantities for three remote sensors is shown in Figure 2 number of images took place in 1999 due to the launch of Landsat 7. The year 2012 witnessed a dramatic drop in the number of images for Landsat 5 ceased to function, followed by a sharp growth in 2013 caused by the launch of Landsat 8.
In addition, for each image, we removed the cloud, cloud shadow, and snow/ice according to the pixel quality attributes of the Landsat surface reflectance product, which were generated from the Fmask algorithm [38]. The preprocessed Landsat time series was used to extract the surface water extent in the GBA.

Auxiliary Data
The Joint Research Centre of European Commission (JRC) Monthly Water History (v1.2) data (https://esdac.jrc.ec.europa.eu/, accessed on 28 June 2021), created by Pickens et al. [3], were used as ancillary data for comparison with the surface water extent extracted by MIWDR [26] in the GBA.
The impervious land surface data were used in this study to represent the extent of urbanization. The data is from the annual maps of global artificial impervious area (GAIA) between 1985 and 2018 by Gong et al. [39,40] (http://data.ess.tsinghua.edu.cn/gaia.html, accessed on 28 June 2021). In addition, the mean overall accuracy of GAIA over multiple years (1985,1990,1995,2000,2005,2010, and 2015) is higher than 90% and the spatial resolution of GAIA is 30 m.
Moreover, annual rainfall data were collected from the Cas Resources and Environment Cloud Platform (https://www.resdc.cn/Default.aspx, accessed on 28 June 2021) to explore the correlation relationship between surface water area and climate change.

Methodology
For clarity and brevity, we established an analytical framework to systematically and quantitively characterize the long-term dynamics of surface water in the GBA for reporting SDG 6.6.1 indicator ( Figure 3). It is divided into two sections: (1) First, we performed data acquisition, preprocessing and surface water extent extraction on the GEE platform to retrieve long time series surface water results for subsequent dynamics analysis. (2) In addition, for each image, we removed the cloud, cloud shadow, and snow/ice according to the pixel quality attributes of the Landsat surface reflectance product, which were generated from the Fmask algorithm [38]. The preprocessed Landsat time series was used to extract the surface water extent in the GBA.

Auxiliary Data
The Joint Research Centre of European Commission (JRC) Monthly Water History (v1.2) data (https://esdac.jrc.ec.europa.eu/, accessed on 28 June 2021), created by Pickens et al. [3], were used as ancillary data for comparison with the surface water extent extracted by MIWDR [26] in the GBA.
The impervious land surface data were used in this study to represent the extent of urbanization. The data is from the annual maps of global artificial impervious area (GAIA) between 1985 and 2018 by Gong et al. [39,40] (http://data.ess.tsinghua.edu.cn/gaia.html, accessed on 28 June 2021). In addition, the mean overall accuracy of GAIA over multiple years (1985,1990,1995,2000,2005,2010, and 2015) is higher than 90% and the spatial resolution of GAIA is 30 m.
Moreover, annual rainfall data were collected from the Cas Resources and Environment Cloud Platform (https://www.resdc.cn/Default.aspx, accessed on 28 June 2021) to explore the correlation relationship between surface water area and climate change.

Methodology
For clarity and brevity, we established an analytical framework to systematically and quantitively characterize the long-term dynamics of surface water in the GBA for reporting SDG 6.6.1 indicator (Figure 3). It is divided into two sections: (1) First, we performed data acquisition, preprocessing and surface water extent extraction on the GEE platform to retrieve long time series surface water results for subsequent dynamics analysis. (2) Second, using the surface water datasets obtained above, we analyzed the changes in surface water from three aspects: 1 The spatial distribution and temporal characteristics of surface water extent were analyzed in the whole GBA and by 11 cities. 2 The annual water inundation frequency time-series for every 30 m×30 m pixel in the GBA was calculated, and dynamic type classification and mapping were carried out based on the temporal variation of the WIF. 3 Combined with impervious land surface data and a land cover transfer matrix, the regions with water losses and gains were identified and analyzed for reporting SDG 6.6.1 indicators. face water from three aspects: ① The spatial distribution and temporal characteristics of surface water extent were analyzed in the whole GBA and by 11 cities. ② The annual water inundation frequency time-series for every 30 m×30 m pixel in the GBA was calculated, and dynamic type classification and mapping were carried out based on the temporal variation of the WIF. ③ Combined with impervious land surface data and a land cover transfer matrix, the regions with water losses and gains were identified and analyzed for reporting SDG 6.6.1 indicators. Figure 3. The framework of this study.

Water Extraction Based on Multiple Index Water Detection Rule
In this study, we employed a simple and fast automatic water extraction method, the multiple index water detection rule (MIWDR), developed by Deng et al. [26], to retrieve surface water time series during 1986-2020 in the GBA for subsequent analysis. The method utilized multiple indices, including NDVI, EVI, MNDWI, and AWEI. Equations (1)-(5) of these indices are as follows:  Figure 3. The framework of this study.

Water Extraction Based on Multiple Index Water Detection Rule
In this study, we employed a simple and fast automatic water extraction method, the multiple index water detection rule (MIWDR), developed by Deng et al. [26], to retrieve surface water time series during 1986-2020 in the GBA for subsequent analysis. The method utilized multiple indices, including NDVI, EVI, MNDWI, and AWEI. Equations (1)-(5) of these indices are as follows: where ρ N IR , ρ Red , ρ Blue , ρ Green , ρ SW IR1 and ρ SW IR2 represent the values of the near infrared band, red band, green band, shortwave infrared band1 and band2 in the Landsat Images. If a pixel meets the following criteria: (AWEI nsh −AWE Ish > −0.1) and (MNDWI > NDVI or MNDWI > EVI), then it is classified as a water body, if not, it is classified as a nonwater body. The MIWDR could be implemented quickly through JavaScript codes based on the GEE cloud computing platform. Compared to other relatively complex machine learning methods, this method has the advantages of simplicity, efficiency, high automation, the capacity to quickly acquire the latest surface water information for analysis.
In addition, we assessed the accuracy of surface water results based on MIWDR method by comparing with the actual visual interpretation results from Landsat true color images. The sample points were created randomly (300 random sample points were selected for each imagery in the study) in each Landsat image tile (10 images in total) via the Create Random Points tool in ArcGIS. These sample points were classified into water and non-water types according to visual interpretation, then a confusion matrix was calculated with the extracted surface water results at the same time. The final accuracy assessment result was presented through overall accuracy and kappa coefficient [21].
Furthermore, we also compared the extracted surface water results with the existing water dataset-the Joint Research Centre of European Commission (JRC) Monthly Water History (v1.2) data in the GBA. The relative error (RE) of surface water extent was calculated to measure the accuracy between the two water datasets. The formula to calculate RE is as follows: where S a respects the surface water area of JRC water dataset, and S e is the surface water area of the extracted water dataset using MIWDR.

Change Trend Analysis Method
Regression analysis is an important method to explore the statistical correlation relationship between different variables and is also suitable for studying the long-term change trend of land surface variables. In this study, the ordinary least square linear regression (OLS-LR) model was used to conduct regression analysis on the time variable x and the surface water area variable y, which could be expressed by the following mathematical formula [41]: In the formula, a is a constant; k is the slope of the linear regression equation; x is the time variable, and ε represents the random error value. This study utilized the least square method to fit the linear relationship between a series of surface water areas and time-series values. The parameter k was used to measure the linear change trend of timeseries data. When k > 0, it represents an upward trend; otherwise, it denotes a downward trend. Moreover, the value of k represents the change speed of surface water extent per year. The k is calculated as follows: In Formulas (2) and (3), x i represents the time value of year i; y i is the surface water area of year i, and N and n represent the total number of observation data points.
Moreover, the goodness of fit test is assessed by using the decisive coefficient R 2 . The closer R 2 is to 1, the closer the whole time-series dataset is to the linear regression line, and the higher the goodness of fit. We applied the OLS-LR over the temporal change analysis in the entire study area as well as 11 cities in the GBA to investigate the surface water extent change trends.

Dynamic Type Classification
Before performing surface water dynamic type classification, we calculated the annual water inundation frequency (WIF) time-series from 2000 to 2020 for each 30 m × 30 m pixel across the study area, considering the Landsat images' quantity over time ( Figure 2). The WIF index directly indicates the frequency of water body occurrence during a certain period of time. The calculation formula of WIF [21] is as follows: where ε i represents the corresponding pixel value of the ith water map, where 1 is water and 0 is non-water; is the total number of observations of water body occurrence, and N represents the total number of valid observations at a specific pixel. Based on the annual WIF time-series of the GBA obtained above, the long time-series dynamic types during 2000-2020 were classified pixel-by-pixel using the mean and range. The thresholds' setting of mean and range was referred to Pickens et al. [42]. The detailed dynamic type classification process for every pixel is shown in Figure 4. Firstly, to reduce annual anomalies, the yearly WIF time-series were smoothed through a three-year average moving window. Secondly, the range and mean of the annual WIF time-series data of each pixel over the period were calculated. The range is the difference between the maximum and minimum value of WIF. In addition, the amount of local extreme values reaching 30% of the range was calculated as the value of WIF change segments. Meanwhile, the permanent and dynamic types of surface water were visualized in RGB color space is visualized according to the WIF change points. The WIF values at the starting position, ending position, and turning position were taken as the values of the red, green and blue channels, respectively. Finally, we achieved surface water dynamic type classification and visualized the spatial distribution.

Identifying Water Losses and Gains under Urbanization
Considering that the GBA is a typical area of rapid urbanization, this study emphasizes the impact of urbanization on the losses and gains of surface water in the study area. Therefore, we only calculated the conversion between surface water bodies and impervious land surfaces through the land cover transition matrix. The land cover transition matrix is the application of the Markov model in land use and land cover, which can quantitatively measure the transformation of different land use types. Furthermore, in order to analyze the spatiotemporal variation in surface water changes by urbanization, we selected seven periods, including 1986, 1990, 1995, 2000, 2005, 2010 and 2015, to calculate the

Identifying Water Losses and Gains under Urbanization
Considering that the GBA is a typical area of rapid urbanization, this study emphasizes the impact of urbanization on the losses and gains of surface water in the study area. Therefore, we only calculated the conversion between surface water bodies and impervious land surfaces through the land cover transition matrix. The land cover transition matrix is the application of the Markov model in land use and land cover, which can quantitatively measure the transformation of different land use types. Furthermore, in order to analyze the spatiotemporal variation in surface water changes by urbanization, we selected seven periods, including 1986,1990,1995,2000,2005,2010 and 2015, to calculate the corresponding surface losses and gains under urbanization.

Accuracy of Extracted Surface Water Bodies in the GBA
A total of 10 images with little cloud coverage (less than 5%) and good quality were selected for accuracy evaluation. These images are evenly distributed every 2-3 years from 1993 to 2020 and cover varied months. As seen from Table 1, the RE (relative errors) of the water area extracted in this study and the JRC water products in the same months are small, ranging from 0.5% to 7.5%, with an average of 4.07%. This indicates that the water area extracted in this study is relatively consistent with the results of JRC products. In addition, the overall accuracy and kappa coefficient results of these 10 images are shown in Table 2. The overall accuracy ranged from 89.95% to 96.89%, with an average value of 93.42%, while the kappa coefficient varied between 0.75 and 0.83, with an average value of 0.794, which indicates the accuracy of the surface water extraction results in this study is relatively high. The spatial distribution of the WIF in the GBA during 1986-2020 is shown in Figure 5, it was divided it into 6 levels: 5-10%, 10-20%, 20-40%, 40-60%, 60-80% and 80-100%. The higher the WIF value is, the more likely the pixel is to be a permanent water body. The   When zooming in some typical regions (a, b, c, d and e), it could be seen that the high-WIF area (>80%) tends to be distributed in the center of rivers and lakes, and the WIF in the marginal region of rivers and lakes shows a decreasing trend from inside to outside (Figure 6e). The spatial distribution characteristics are consistent with the natural and hydrological principles of river and lake dynamics, and it is well known that the impact of human activities on water bodies often occurs around surface water Figure 6a,b). In addition, it should be noted that in the interlaced area of river networks, a substantial number of water bodies are distributed in blocks or fragments (Figure 6c), and their WIF values are mostly less than 50%, which could be attributed to the abundance of paddy fields and pits in the region.

Temporal Changes in Surface Water in the GBA from 1986 to 2020
For the whole GBA, the annual surface water body area calculated under different WIF thresholds from 1986 to 2020 is shown in Figure 7. The permanent surface water bodies (WIF > 75%) did not change much. In contrast, the maximum surface water extent (WIF > 0%) in the GBA demonstrates a significant decreasing tendency after 2001 (R 2 = 0.8566, p < 0.0001), with a loss of 155.41 km 2 per year. Furthermore, the seasonal surface water extent (35% < WIF < 75%) fluctuates slightly over the last 35 years. It also could be seen that the lower the WIF threshold is, the higher the loss rate after 2001.

Temporal Changes in Surface Water in the GBA from 1986 to 2020
For the whole GBA, the annual surface water body area calculated under different WIF thresholds from 1986 to 2020 is shown in Figure 7. The permanent surface water bodies (WIF > 75%) did not change much. In contrast, the maximum surface water extent (WIF > 0%) in the GBA demonstrates a significant decreasing tendency after 2001 (R 2 = 0.8566, p < 0.0001), with a loss of 155.41 km 2 per year. Furthermore, the seasonal surface water extent (35% < WIF < 75%) fluctuates slightly over the last 35 years. It also could be seen that the lower the WIF threshold is, the higher the loss rate after 2001. The yearly surface body areas of eleven individual cities in the GBA during 1986-2020 showed divergent trends over time (Figure 8). For better separation, we classified these cities into three groups to illustrate their change characteristics according to the range of surface water area. Furthermore, the linear trend lines were added to the line charts corresponding to each city. Among the cities, the surface water areas in Guangzhou, Foshan, Huizhou, Zhuhai, Shenzhen, Hong Kong and Macao show an evident downward trend. Especially for Macao, the linear declining trend in Macao is very significant, with R 2 of 0.96, but its loss rate is fairly low, with about 0.426 km 2 per year. It is noteworthy that there are no cities with a significant upward trend here, which reflects that the majority of cities are experiencing a decrease in surface water areas to varying degrees. Therefore, for cities with a rapid downward trend, more attention should be given to surface water protection and recovery. The yearly surface body areas of eleven individual cities in the GBA during 1986-2020 showed divergent trends over time ( Figure 8). For better separation, we classified these cities into three groups to illustrate their change characteristics according to the range of surface water area. Furthermore, the linear trend lines were added to the line charts corresponding to each city. Among the cities, the surface water areas in Guangzhou, Foshan, Huizhou, Zhuhai, Shenzhen, Hong Kong and Macao show an evident downward trend. Especially for Macao, the linear declining trend in Macao is very significant, with R 2 of 0.96, but its loss rate is fairly low, with about 0.426 km 2 per year. It is noteworthy that there are no cities with a significant upward trend here, which reflects that the majority of cities are experiencing a decrease in surface water areas to varying degrees. Therefore, for cities with a rapid downward trend, more attention should be given to surface water protection and recovery.

Dynamic Type Mapping and Analysis in the GBA during 2000-2020
The dynamic type map for surface water variations during 2000-2020 in the GBA is illustrated in Figure 9, in which the color represents varied dynamic types (seen in Figure  4). As seen in some typical enlarged areas of dynamic types, the edge of the lake in subregion A is orange, and the central area of the lake is white, which indicates the surface

Dynamic Type Mapping and Analysis in the GBA during 2000-2020
The dynamic type map for surface water variations during 2000-2020 in the GBA is illustrated in Figure 9, in which the color represents varied dynamic types (seen in Figure 4). As seen in some typical enlarged areas of dynamic types, the edge of the lake in subregion A is orange, and the central area of the lake is white, which indicates the surface water around the lake is losing. In contrast, the blue edge of the lake in subregion B presents an expansion trend. Some patches experiencing a dry period (with WIF growing down first and up later) in paddy fields in the subregion C are pink, while the green region of the meandering river streamway (in the subregion D) denotes that it witnessed a wet period in history. The lakes with light gray edges in subregion E experience seasonally stable changes, which could be attributed to the regular hydrological phenology. The paddy fields in subregion F exhibit signs of high-frequency changes, which may be caused by intensive human activities. The statistics of dynamic type classification results are presented in Table A1 in detail and visualized in Figure 10. All water-related areas account for 6.48% (3594.01 km 2 ) of the total area in the GBA, of which 57.46% was dynamic types during 2000-2020. This suggests that most surface water bodies change over the last 20 years instead of remaining unchanged. In addition, the gain, loss, dry and wet periods account for approximately 4-7% of all areas with water, while stable seasonal and high-frequency periods account for approximately 16-18% of all areas with water. The difference in areas indicates that the surface water area experiencing multi-direction changes is greater than that of single-direction changes. In reality, to support the achievement of SDG 6.6.1, areas of loss should be given more attention from local governments and residents. The statistics of dynamic type classification results are presented in Table A1 in detail and visualized in Figure 10. All water-related areas account for 6.48% (3594.01 km 2 ) of the total area in the GBA, of which 57.46% was dynamic types during 2000-2020. This suggests that most surface water bodies change over the last 20 years instead of remaining unchanged. In addition, the gain, loss, dry and wet periods account for approximately 4-7% of all areas with water, while stable seasonal and high-frequency periods account for approximately 16-18% of all areas with water. The difference in areas indicates that the surface water area experiencing multi-direction changes is greater than that of singledirection changes. In reality, to support the achievement of SDG 6.6.1, areas of loss should be given more attention from local governments and residents.
unchanged. In addition, the gain, loss, dry and wet periods account for approximately 4-7% of all areas with water, while stable seasonal and high-frequency periods account for approximately 16-18% of all areas with water. The difference in areas indicates that the surface water area experiencing multi-direction changes is greater than that of single-direction changes. In reality, to support the achievement of SDG 6.6.1, areas of loss should be given more attention from local governments and residents.

Surface Water Losses and Gains under Urbanization in the GBA
The overall spatial distribution of surface water losses and gains is shown in Figure 11A,B. In the two enlarged regions (Figure 11f,g) of water losses, it was found that many cities and towns occupy water bodies in coastal areas and the reclamation of coastal urban zones. In the zoomed region (Figure 11h,i) of water gains, some built-up areas are transformed into water due to the construction of artificial reservoirs and lake parks, like Gongming reservoir and Wenhan Lake Park. In general, the area of water bodies occupied by urbanization is far greater than that restored to water bodies.
Through the land cover confusion transfer matrix, the area of water bodies transformed into impervious surfaces and the area of built-up areas transformed into water bodies in different years were calculated( Figure 12). The surface water losses area shows a significant declining trend from 211.29 km 2 in 1986 to 2.92 km 2 in 2015 (Figure 12a), and the surface water restoration area shows a slightly growing trend (Figure 12b), which indicates that more water bodies have been protected in recent years and that more attention has been given to water restoration.
The overall spatial distribution of surface water losses and gains is shown in Figure  11A,B. In the two enlarged regions (Figure 11f,g) of water losses, it was found that many cities and towns occupy water bodies in coastal areas and the reclamation of coastal urban zones. In the zoomed region (Figure 11h,i) of water gains, some built-up areas are transformed into water due to the construction of artificial reservoirs and lake parks, like Gongming reservoir and Wenhan Lake Park. In general, the area of water bodies occupied by urbanization is far greater than that restored to water bodies. Through the land cover confusion transfer matrix, the area of water bodies transformed into impervious surfaces and the area of built-up areas transformed into water bodies in different years were calculated( Figure 12). The surface water losses area shows a significant declining trend from 211.29 km 2 in 1986 to 2.92 km 2 in 2015 (Figure 12(a)), and the surface water restoration area shows a slightly growing trend (Figure 12(b)), which indicates that more water bodies have been protected in recent years and that more attention has been given to water restoration.

Possible Causes for Surface Water Extent Changes
Climate change is an important factor that may affect surface water area. Precipitation is a key component of climatic factors, which has a complicated relationship with surface water extent changes. We briefly analyzed the relationship between precipitation

Possible Causes for Surface Water Extent Changes
Climate change is an important factor that may affect surface water area. Precipitation is a key component of climatic factors, which has a complicated relationship with surface water extent changes. We briefly analyzed the relationship between precipitation and maximum surface water extent by comparing the changes during 1986-2019. We also calculated the linear correlation coefficient between them ( Figure 13). The change trend of precipitation and surface water extent was inconsistent and the linear correlation between them is not significant (p < 0.05), but the change trend was relatively consistent over some short-term periods, which may be caused by extreme precipitation events.

Possible Causes for Surface Water Extent Changes
Climate change is an important factor that may affect surface water area. Precipitation is a key component of climatic factors, which has a complicated relationship with surface water extent changes. We briefly analyzed the relationship between precipitation and maximum surface water extent by comparing the changes during 1986-2019. We also calculated the linear correlation coefficient between them ( Figure 13). The change trend of precipitation and surface water extent was inconsistent and the linear correlation between them is not significant (p < 0.05), but the change trend was relatively consistent over some short-term periods, which may be caused by extreme precipitation events.

Comparison with Other Studies
Comparing the surface water area derived in this study with the result of other researchers, it is found that the results are highly consistent. According to the study of Zhao Y L, et al. [43], the surface water area of the GBA in 2015 was 5266.02 km 2 , which was obtained based on remote sensing data used in the Second National Land Survey, China, in 2015 and checked for correctness. While the maximum surface water area of the GBA Figure 13. Correlation analysis between mean annual precipitation and maximum surface water area during 1986-2019 in the GBA.

Comparison with Other Studies
Comparing the surface water area derived in this study with the result of other researchers, it is found that the results are highly consistent. According to the study of Zhao Y L, et al. [43], the surface water area of the GBA in 2015 was 5266.02 km 2 , which was obtained based on remote sensing data used in the Second National Land Survey, China, in 2015 and checked for correctness. While the maximum surface water area of the GBA in 2015 was 5318.80 km 2 in this study, the difference between the two results was 52.78 km 2 , which further verified the reliability of the surface water results in this study. However, we are not able to compare the temporal variations due to a lack of relevant studies in the GBA.
In addition, the surface water monitoring and analysis method and framework of this study features by efficiency of the GEE platform and serving the goal of SDG 6.6.1. The most important point is that we analyzed the dynamic types of surface water to characterize and visualize the spatiotemporal variations for reporting SDG 6.6.1 in the study area. In comparison with previous studies related to surface water extraction and SDG 6.6.1 [2,21], this study's analytical framework is more refined and comprehensive.

Implications for SDG 6.6.1 Reporting and Policy-Making
Since the SDG 6.6.1 indicator aims to track variations in different types of water-related ecosystems, informing decision-makers of the extent changes of water-related ecosystems over time, the workflow, framework and analysis of this study are highly consistent with these targets, which might be taken as references for other similar studies related to SDG 6.6.1 reporting and policy-making of regional government for sustainable development.
Additionally, from the surface water extent spatiotemporal variations derived in this study, it is evident that the surface water area in the GBA shows a slow downward trend in recent 35 years and the surface water area keeping losing occupies 7.34% (263.75 km 2 ) of all water-related area in the GBA. Therefore, more attention should be paid to regions experiencing decreasing trends to pre-empt further loss in surface water resources.

Limitations and Future Improvements
There are still some limitations and uncertainties in this study.
(1) The Landsat imagery used in this study is easily affected by clouds, and the GBA is a typical coastal region with frequent cloud coverage. This issue limits the highly frequent monitoring of surface water when it comes to extreme disasters (e.g., extreme precipitation events) in the GBA. (2) In regard to surface water verification, we selected ten images for accuracy assessment, the number of images is relatively small and the evaluation of images with a medium or high cloud coverage was not taken into consideration, which could cause some uncertainties to the surface water results. (3) The exploration of driving factors of long-term surface water variations is still not enough. We only analyzed the impact of land cover change and precipitation on the surface water changes. Due to low accessibility to more meteorological and hydrological in the GBA and the complexity of mutual influence mechanism between human activities and climate changes, we did not use more datasets to make a comprehensive analysis of the driving forces.
In the future research, we will optimize the data preprocessing procedures to reduce the effects of clouds on the optical remote sensing imagery and introduce other multi-source remote sensing data, such as Sentinel 1/2, HJ, GF and other data with high resolution and quality to better monitor the spatiotemporal dynamics of surface water in the study area. Additionally, we need to apply multi-factor analysis to investigate the driving forces of surface water changes in the GBA more comprehensively. Moreover, we are supposed to integrate.

Conclusions
In this study, in order to monitor SDG 6.6.1 indicator, we utilized long-term Landsat images and the MIWDR method to extract surface water extent of the GBA during 1986-2020 on the GEE platform. We characterized the spatial and temporal variations in surface water in the GBA and classified dynamic types of surface water by pixels to explore the features of surface water dynamics in detail. In addition, the surface water losses and gains caused by urbanization were identified and analyzed. The main conclusions are as follows: (1) From 1986 to 2020, the average minimal and maximal surface water extents of the GBA were 2017.62 km 2  year, while the surface water restored from impervious land surfaces denoted a slightly growing trend. It appears that surface water extent is under less pressure induced by urbanization in recent years and has more chances to be restored and protected by human activities.
In addition, we also assessed the water extraction accuracy and possible climate factors for surface water changes. However, due to the uncertainties and limitations that existed in our study, we need to further explore other factors affecting the surface water extent of GBA in following research.
In this study, the long-term surface water time-series datasets and spatiotemporal analysis results could be provided for reference by local governments or related ecological environment protection agencies to assess regional procedures toward achieving SDG 6.6.1 as well as to inform environmental protection policy-makers of the key areas for water conservation and restoration in the GBA. Furthermore, the methodology, analytical framework and remote sensing big data resources available on the GEE could be easily applied in other similar regions worldwide for supporting the SDG 6.6.1 indicator and promoting regional water resources protection as well as management.  Acknowledgments: We want to express our respect and gratitude to Chunlin Wang, Xiaoya Wang, Zhuo Li, Yuhao Pan and Xinhao Pan for their help. We also would like to thank the reviewers for their valuable comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.