Mapping of the Spatial Scope and Water Quality of Surface Water Based on the Google Earth Engine Cloud Platform and Landsat Time Series

Abstract

Surface water is an important parameter for water resource management and terrestrial water circulation research that is closely related to human production and livelihood. With the rapid development of remote sensing technology and cloud computing platforms, the use of remote sensing technology for large-scale and long-term surface water monitoring and investigation has become a research trend. Based on the Google Earth Engine (GEE) cloud platform and Landsat series satellite data, in this study, the Emergency Geomatics Service (EGS) operational surface water mapping algorithm and water index masking were utilized to extract the spatial scope of the water body. The validated models of the Secchi disk depth (SDD), chlorophyll-a (Chl-a) and suspended solids (SS) concentration were applied to water quality parameter inversion and water quality evaluation. Surface water extent extraction and water quality maps were created to analyze the spatial distribution of the water body and the spatial–temporal evolution characteristics of the water quality parameters. A verification experiment was carried out with the surface water in Zhejiang Province as the research object. The results show that the surface water in the study area from 1990 to 2022 could be accurately extracted. The kappa coefficients were all greater than 0.90, and the overall accuracies of the extractions were greater than 95.31%. From 1990 to 2022, the total surface water area in Zhejiang Province initially decreased and then increased. The minimum water area of 2027.49 km² occurred in 2005, and the maximum water area of 2614.96 km² occurred in 2020, with an annual average variation of 193.92 km². Since 2015, the proportion of high SS and Chl-a concentrations, and low SDD water bodies in Zhejiang Province have decreased, and the proportion with better water quality has increased significantly. The spatial distribution map of the surface water and the inversion results of the water quality parameters obtained in this study provide a valuable reference and guidance for regional water resource management, disaster monitoring and early warning, environmental protection, and aquaculture.

1. Introduction

Surface water refers to the water that exists on the surface of the Earth’s crust and is exposed to the atmosphere, including rivers, lakes, reservoirs, and ponds. It plays a crucial role in the balance between ecosystems and socio-economic development and provides essential resources for the survival of living organisms and the development of human society [1,2,3]. Surface water resources have a variety of ecosystem service values such as water storage, flood control, aquaculture, power generation, and transportation. With the gradual shortage of water resources, the social function of surface water is becoming increasingly significant, and the economic importance of surface water is continuously increasing [4,5]. The area of surface water resources has shrunk, and the pollution of the water environment has substantial negative impacts on the ecological environment and human development [6,7,8]. Therefore, understanding the spatial extent and time series features of surface water resources and investigating the variations in water quality parameters are significant for ecological protection, social production, and life development.

With the progress of remote sensing technology and the rapid development of cloud computing technology, the traditional hardware environment dependence problem has been solved, which makes large-scale and long-term water body monitoring more efficient [9]. Based on the characteristics of water spectral curves, scholars have proposed a variety of water information extraction methods, such as the single band threshold-based method, the spectrum photometric-based method, and the water index-based method. The single-band threshold-based method uses the characteristics of water’s strong absorption and low spectral reflectance in the near-infrared and shortwave infrared bands to distinguish water from other ground objects by setting a threshold [10]. This method has the benefits of simplicity of operation and easy implementation. However, it does not make full use of multi-band information and is often affected by other ground objects with similar spectral characteristics. The spectrum photometric-based method identifies the spectral characteristics of water information by analyzing the spectral characteristic curves of ground objects in the original band or transformation space of remote sensing images. Then, it uses logical judgment rules to extract water information [11]. The spectrum photometric-based method has a high precision and wide application and has good adaptability to the complex distribution of ground objects. Nevertheless, the model construction process is cumbersome and the computational complexity is high [12]. The most commonly used water index-based method is based on the ratio operation between characteristic bands to enhance the contrast between water and other features, which is used to identify water information. Mcfeters (1996) first proposed the normalized difference water index (NDWI), using the ratio of band 2 (green band) to band 4 (near-infrared band) of Landsat thematic imager (TM) to recognize the water body. Xu (2006) used Landsat TM bands 2 (green band) and 5 (mid-infrared band) to evaluate the modified NDWI (MNDWI) [13,14]. Feyisa et al. developed a new index (the automated water extraction index, AWEI) using the five spectral bands of the Landsat 5 TM to increase the contrast between water and other dark surfaces and to maximize the separability of water and non-water pixels through band difference solution, addition, and the application of different coefficients [15]. With the development of computer technology, such as supervised classification, decision tree algorithms, the random forest method, and other machine learning algorithms have also been widely used in remote sensing water information extraction. Tulbure et al. used the random forest model to study the scope and dynamic changes of surface water and floods (1986–2011) in the Murray-Darling Basin (MDB). In addition to band reflectance, the explanatory variables of the model also include the brightness temperature, normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), digital elevation model (DEM) derivative variables, and multiple water body indexes [16]. Wang et al. proposed a multi-dimensional and densely connected new deep convolutional neural network (DenseNet) for identifying water bodies from high-resolution multi-spectral images. Its extraction accuracy is higher than those of the water body index method and classical convolutional neural networks [17].

Water quality parameters refer to indicators that reflect the specific substances in or physical properties of the water, such as the dissolved oxygen content (DOC), pH value, turbidity, total nitrogen content, and total phosphorus content. Water quality parameters provide data support for the evaluation of water quality conditions, and shift the evaluation criteria from subjective judgment to quantitative description, which is more objective, scientific, and reliable. For water color remote sensing, any water component parameter with significant spectral characteristics or optical properties, namely, photoactive compounds, can be quantitatively inverted via remote sensing, such as the chlorophyll content, suspended solids, and DOC. The inversion methods mainly include empirical/semi-empirical methods and analytical/semi-analytical methods [18]. The empirical/semi-empirical model mainly realizes the inversion of water quality parameters based on the statistical relationships between the water quality parameters and remote sensing parameters. The semi-empirical method adds the theoretical analysis of the relationships between remote sensing parameters and water quality parameters to the statistical relationship. The core of the analytical/semi-analytical method is the bio-optical model. Thus far, numerous experts and scholars have constructed models using different water quality parameters based on Landsat data. Based on Landsat data and actual measured water quality data, Brivio et al. used regression analysis to establish an inversion model of the chlorophyll-a concentration in Lake Garda [19]. Sharaf El Din et al. extracted the principal components to develop algorithms for turbidity and total suspended solids (TSS) based on Landsat-8 surface reflectance data [20]. According to Landsat operational land imager (OLI) data, Song et al. established a remote sensing model for estimating transparency using the ratio of the red to blue bands and achieved the inversion of the global water transparency [21].

The main purpose of this study is to extract the surface water coverage based on the Google Earth Engine (GEE) cloud platform and long-term Landsat remote sensing images and to retrieve water quality parameters, including chlorophyll-a (Chl-a) and suspended solids (SS) concentration, and transparency, which are characterized by the Secchi disk depth (SDD). The spatial extent of the surface water and the spatial–temporal distribution of the water quality in Zhejiang Province from 1990 to 2022 were mapped, and the characteristics of the spatial–temporal variations in the water area and water quality parameters were analyzed. The results of this study provide data and technical support for the scientific management of regional water resources.

2. Study Area and Datasets

2.1. Study Area

Zhejiang Province, China, is located in the south wing of the Yangtze River Delta along the southeast coast of China (27°02′–31°11′N, 118°01′–123°10′E). It has a land area of 105,500 km². The terrain of Zhejiang Province is diverse, mainly including mountains and hills. The terrain elevation decreases from southwest to northeast in a stair-step manner. Zhejiang Province is located in the subtropical monsoon zone and has a mild and humid climate, abundant water resources, and a well-developed water system. The geographical location and topographic map of the study area is shown in Figure 1.

2.2. Datasets

Landsat is an earth observation satellite mission jointly planned by National Aeronautics and Space Administration (NASA) and the United States Geological Survey (USGS). The project provides valuable data for a wide range of applications, including land use planning, resource management, and environmental monitoring. All Landsat images are publicly available through the USGS Earth Explorer. The 30 m resolution remote sensing images (including surface reflectance data and top of atmosphere data) of Zhejiang Province from 1990 to 2022 originated from the remote sensing data of the Landsat 5, 7, and 8 satellites. The European Commission Joint Research Centre (JRC) provides a rich, global long-term surface water data set, which provides data support for the investigation and research of water resources [22]. The spatial and temporal distribution information of global surface water from 1984 to 2021, as well as statistical data about the scope and variation, such as the occurrence (the frequency with which water was present), was obtained from the JRC global surface water dataset. The nearly global elevation data were obtained from the Shuttle Radar Topography Mission (SRTM) global digital elevation model released by NASA. Due to the stability of the elevation data and the span of the study period, in this study, we only used the digital elevation data for 2000 instead of data for all of the research years. Landsat Global Inland Water, a global inland water body dataset based on Landsat satellite data provided by GLCF (Global Land Cover Facility), includes the changes in inland water bodies around the world over the last 35 years. The datasets we used in this study are also accessible on the GEE platform [23]. The first 1 m resolution national-scale land-cover map of China, SinoLC-1, was established using a deep learning-based framework and open-access data including global land-cover (GLC) products, open street map (OSM), and Google Earth imagery [24]. The SinoLC-1 was created by the Low-to-High Network (L2HNet), which can be found at: https://www.sciencedirect.com/science/article/abs/pii/S0924271622002180?via%3Dihub (accessed on 15 March 2023). The datasets used in this study is shown in

Table 1. Landsat time series and GEE datasets.

Landsat Time Series		Datasets
Years	Satellite Sensor	Years	Data	Resolution
1990–1995 2000–2010 2015–2022	Landsat5 TM Landsat7 ETM+ Landsat8 OLI	1984–2021	JRC Global Surface Water Mapping Layers, v1.4	30 m
		2000	NASA SRTM DEM
		2000	GLCF: Landsat Global Inland Water
		2021	SinoLC-1	1 m

.

3. Methods

Based on the GEE cloud platform and long-term Landsat remote sensing data, in this study, the Emergency Geomatics Service (EGS) operational surface water mapping algorithm was used to extract the surface water coverage area, and the verified model was utilized to retrieve the Chl-a, SS, and SDD. Water body distribution and water quality maps were created to analyze the temporal and spatial variation characteristics of the water body area and water quality parameters. The technical workflow of this study is shown in Figure 2.

3.1. Data Pre-Processing

In this study, the top of atmosphere (TOA) data and surface reflectance (SR) data from Landsat data, which had been geometrically corrected, were obtained from the (GEE) cloud platform. The existence of clouds has a great impact on remote sensing information extraction and target recognition. To ensure the quality of the remote sensing images, TOA data with less than 20% cloud coverage were selected when filtering data, and the quality assessment band was used for the cloud removal processing to generate long-term Landsat image data from 1990 to 2022 with an interval of 5 years using the median synthesis method. Based on the Zhejiang border shape file, the remote sensing images were cropped and the TOA data and SR data covering the study area were retained.

3.2. Water Spatial Range Mapping Based on EGS Operational Surface Water Mapping Algorithm

The EGS operational surface water mapping algorithm is a direct and fully automatic surface water extent extraction method. This algorithm uses historical inundation frequency products to mark land and water areas and sample the feature information; then, it implements land and water classification using the random forest algorithm to realize the extraction of the surface water [21]. Random forest is an ensemble algorithm, which belongs to the bagging type. It takes the decision tree as the basic unit and forms a random forest by integrating a large number of decision trees. Each time, it samples with replacements from the training set to form a new training set, and uses the new training set to train sub-models. The final results are obtained by voting or averaging, so that the results of the overall model have higher accuracy and generalization performance.

The JRC Global Surface Water Mapping Layers v1.4 dataset obtained from the GEE cloud platform contains the location and temporal distribution information about global surface water from 1984 to 2021 and includes statistical data on the scope of and changes in these water surfaces [25]. For the purpose of ensuring that all of the water sampling points are water samples, the algorithm generates permanent water data (water occurrence frequency of >80%) based on JRC global surface water data for masking. It extracts six bands of information (blue, green, red, near-infrared, and shortwave infrared 1 and 2) from the TOA data using the permanent water mask, combined with the permanent land (water occurrence frequency of 0%) as training samples. The number of samples in each category was greater than 600, with a ratio of 7:3 for the training and validation datasets. In the random forest algorithm, the number of decision trees was set to 100 and the separation of water and land was realized.

When only using band attributes for random forest training to extract the water bodies, the characteristics of the same spectrum of remote sensing images have a certain impact on the extraction results. For example, in this study, some mountain shadows were misjudged as water bodies. In order to further improve the classification accuracy and eliminate the influence of misclassification and singular values, the MNDWI and NDWI were used as masks to enhance the water body information. Slope calculation was conducted based on NASA’s SRTM DEM (30 m), and the data with slopes of <20° were selected for masking. According to the high vegetation coverage characteristic in the shadowed part of the mountains, which indicates that the NDVI value was also high, we used the NDVI as a mask to remove the misclassification caused by the shadows of the mountains. The comparison of surface water extent extraction efficiency using RF method is shown in Figure 3.

3.3. Accuracy Evaluation of Surface Water Extent Extraction Results

In this study, the accuracy evaluation of surface water extent extraction results is divided into two parts. The first is the accuracy evaluation of the random forest training model, which is measured by two indexes of kappa and overall accuracy (OA) [26]. The Kappa coefficient is an indicator for consistency tests and can also be used to measure the effect of classification. Because of the classification problem, consistency is based on whether the model prediction results are consistent with the actual classification results. The calculation of the kappa coefficient is based on the confusion matrix, and the value is between −1 and 1.

The kappa coefficient calculation formula based on the confusion matrix is as follows:

$$Kappa=\frac{p_o-p_e}{1-p_e}$$

where

$$p_0=\frac{Sum of diagonal elements}{ \mathrm{S}um of the whole matrix elements},$$

$${\mathrm{p}}_e=\frac{{\displaystyle \sum _i^{}Sum of the elements in row i\times Sum ofthe elements in column i }}{({\displaystyle \sum _{}^{}All elements of the matrix{)}^2}}$$

Overall accuracy: the sum of correctly classified pixels divided by the total number of pixels is used to describe the overall classification accuracy.

$$OA=\frac{N_{\mathrm{a}}}{N_{\mathrm{t}}}$$

where $N_{\mathrm{a}}$ is the number of correctly classified pixels, and $N_t$ is the total number of pixels.

The greater the kappa and OA index values of the results obtained by statistical calculation, the better the model classification results and the higher the accuracy. In addition, the outcomes of this study are compared with the other two surface water data products and remote sensing source images in the same region. The accuracy and error of surface water extent extraction can be more intuitively and clearly shown through the display of the effect diagrams in small areas.

3.4. Inversion and Mapping of Water Quality Parameters

Water resources are an indispensable natural resource for humans and are a fundamental condition for human survival and development. However, with the development of society, the pollution of water resources and the environment has become increasingly serious. As a basic component of water pollution control and treatment, water quality detection has become vitally significant. Water conditions can be evaluated using various water quality parameters, including the water transparency, suspended solids concentration, chlorophyll concentration, and dissolved organic matter. Therefore, the inversion of water quality parameters is of great significance for water resource management, water environment protection, and water pollution prevention and control. Regarding water quality parameter inversion methods, the empirical/semi-empirical algorithm is generally more accurate and convenient to use, but it lacks a certain physical basis. The physical basis of the analysis model is powerful; nonetheless, the overly complex physical mechanism leads to the demand for a large amount of field data and cumbersome operation processing in the inversion of the model. In this study, a large number of previous studies were selected based on spatial–temporal conditions and some actual site data. Based on the water mask after the land–water separation, the Landsat SR data were masked and combined with the screened models to retrieve three water quality parameters in Zhejiang Province, including the Chl-a, SS, and SDD. The models of Chl-a, SS and SDD are shown in

Table 2. Information about the water quality parameter inversion models.

Parameters	Models	References
Chl-a	$Y_{Chl-a}=124.3\times R_{red}{}^2+15.28\times R_{red}+0.914$	[27]
SDD	$Y_{SDD}=6715.539\times {({\mathrm{R}}_{blue}+{\mathrm{R}}_{green})}^2-2794.743\times ({\mathrm{R}}_{blue}+{\mathrm{R}}_{green})+336.465$	[28]
SS	${\ln }_{SS}=14.656\times ({\mathrm{R}}_{green}+{\mathrm{R}}_{red})/({\mathrm{R}}_{green}/{\mathrm{R}}_{red})+1.661$	[29]

.

3.5. Analysis of the Spatial Distribution of the Water Bodies and the Evolution of the Temporal-Spatial Patterns of the Water Quality Parameters

Based on the spatial mapping results of the surface water, statistical analysis of the area in the entire region and in different districts and cities in Zhejiang Province in different years was conducted, and the temporal variations in the spatial distribution of the water bodies were obtained. In order to characterize the degree of the dynamic changes in the surface water area in each region, the coefficient of variation (CV) of each district was calculated separately.

$$C{V}_{(i)}=\frac{{\mathrm{\mu }}_i}{\overline{x_i}},$$

(1)

where ${\mathrm{\mu }}_i$ is the standard deviation, and $\overline{x_i}$ is the average value. If the CV value is high, there is a large difference between the individual data and the overall average value (or central tendency), and the data fluctuation is more intense; in other words, the water area in the region changes greatly.

Grade division was performed on the inversion results of the water quality parameters according to the different concentration ranges. The proportion of water bodies with different concentrations by year was compiled to obtain the spatial–temporal variation trends of the different concentration water bodies under different water quality parameters. Through calculations and drawing the Sankey diagram for the different concentrations in the same research phases, the transformation of the water bodies with different concentrations was vividly depicted.

4. Results and Analysis

4.1. Spatial Distribution Mapping of Surface Water

The spatial distribution characteristics of the surface water in each region and the extraction efficiency of the water body details are shown in Figure 4. The integrity and continuity of the surface water extent extraction were good, and the main river system characteristics in Zhejiang Province were well preserved. The overall shape and edge characteristics of Qiandao Lake (Figure 4D), the largest lake in Zhejiang Province, were effectively extracted. Lakes such as Hangzhou West Lake (Figure 4B) and Ningbo Dongqian Lake (Figure 4F) were accurately and completely extracted. Moreover, the small island in the lake was accurately distinguished and retained. The lakes, ponds, and crisscrossing rivers in the Hangzhou–Jiaxing–Huzhou Plain in northern Zhejiang Province (Figure 4A) and the scattered small water bodies in the Jinqu Basin (Figure 4E) were successfully distinguished and extracted. It was found that the algorithm utilized in this study has the ability to accurately identify the river network within the city, and it retained the continuity characteristics of the river while extracting the water body (Figure 5).

The Kappa coefficients for each year were greater than 0.90, and the OA values were greater than 0.95 (

Table 3. The performance of the surface water extent extraction accuracy for each year.

	1990	1995	2000	2005	2010	2015	2020	2022
OA (%)	0.9531	0.9589	0.9677	0.9795	0.9736	0.9736	0.9648	0.9735
Kappa	0.9061	0.9179	0.9354	0.9589	0.9471	0.9472	0.9295	0.9470

), indicating that the water information extraction results are credible. Furthermore, the surface water extent extraction results were compared with the effect pictures of GLCF, SinLC-1 products, and remote sensing source images in the same small area. In contrast with GLCF products, the products in this study have better effects on river continuity, integrity, and water bodies shape feature extraction. Despite SinLC-1 product with 1 m resolution, they also have high consistency.

Based on the extraction of the water body information in each year, the surface water area of Zhejiang Province was counted with eleven districts as units (

Table 4. Statistics of surface water area in 11 districts of Zhejiang Province.

Area (km2)	1990	1995	2000	2005	2010	2015	2020	2022	$\overline{{\mathit{x}}_{\mathit{i}}}$
Hangzhou	897.96	872.13	828.72	754.06	825.89	873.74	854.20	834.56	842.65
Ningbo	168.02	150.65	162.21	126.08	159.87	172.30	184.95	187.71	163.97
Wenzhou	167.39	150.48	173.69	150.66	159.36	162.66	158.49	160.63	160.42
Jiaxing	167.42	145.16	148.27	113.37	127.42	143.11	150.33	153.58	143.58
Huzhou	268.88	228.04	254.23	203.76	243.47	298.63	317.94	304.89	264.98
Shaoxing	216.26	183.96	253.02	209.88	220.45	252.27	252.35	230.56	227.34
Jinhua	169.91	159.81	178.10	166.71	177.69	219.24	210.76	187.20	183.67
Quzhou	102.84	99.88	111.96	80.73	124.89	154.86	160.41	146.22	122.72
Zhoushan	7.94	5.82	7.09	5.19	4.81	6.78	7.01	7.35	6.49
Taizhou	156.68	139.55	146.02	134.50	149.47	155.42	162.86	163.83	151.04
Lishui	101.13	88.64	106.06	82.56	145.37	171.30	155.65	155.71	125.80
Total area	2424.41	2224.12	2369.3	2027.49	2338.70	2610.30	2614.96	2532.24	2392.69

).

As can be seen from Table 4 and Figure 6, the average surface water coverage area in Zhejiang Province from 1990 to 2022 was 2392.69 km². The total surface water area exhibited a decreasing trend from 1990 to 2005 and reached the minimum value of 2027.49 km² in 2005. Subsequently, the surface water coverage area continuously increased, and the maximum value of 2614.96 km² was reached in 2020. The average annual change in the area was 193.92 km². According to the statistics of the surface water data for eleven districts and cities in Zhejiang Province, the interannual variation in the surface water area was roughly the same as that of the total surface water area in Zhejiang, but there were some differences in individual cities. The average annual surface water area of Hangzhou is 842.65 km², far exceeding that of other districts, while Zhoushan only has 6.49 km², which is only 0.77% of Hangzhou. The surface water area in other districts is between 100–300 km².

4.2. Inversion of Water Quality Parameters and Water Quality Map

4.2.1. Chlorophyll-a Concentration

A Landsat SR median synthetic image for the study year was acquired through the GEE cloud platform; then, the Chl-a concentration was inverted using the validated model. Finally, the spatial and temporal distributions of the Chl-a concentration were mapped using the water mask data (Figure 7).

Based on the data presented in Figure 7, the Chl-a concentrations of the water bodies in Zhejiang Province exhibited a decreasing trend from 1990 to 2022. The distribution pattern indicates that the concentration was lower in the mountainous regions, followed by the hilly and basin areas in the central region, and the highest concentration occurred in the coastal areas. The locations with rapidly increasing Chl-a concentrations strongly coincided with the distributions of several major cities.

To grade the results of the Chl-a concentration inversion, the outliers in the first and last 1% of the dataset were removed, and the water bodies were divided into five different gradations: low (0–1.5 µg/L), medium–low (1.5–3.0 µg/L), medium (3.0–4.5 µg/L), medium–high (4.5–6.0 µg/L), high (6.0–7.5 µg/L). The Sankey diagram was created to visualize the transition between the different grades of the Chl-a concentration (Figure 8).

As can be seen from Figure 8, the percentage of low Chl-a concentrations increased annually, and the proportion of high concentrations increased until 2010 and then rapidly decreased. However, the percentage of low Chl-a concentrations did not exceed 10% in all of the periods studied. The transformation between the low and medium–low Chl-a concentration categories was much stronger than in the high-value region. The concentration variations were mainly between the low, medium–low, and medium categories. A large number of medium–low Chl-a concentration water bodies shifted to medium Chl-a concentration water bodies in the years when the Chl-a concentrations increased significantly, while the medium–low Chl-a concentration water bodies shifted to low Chl-a concentration water bodies in years when they decreased dramatically. A large proportion of the medium–high concentrations changed, but throughout the time series, they remained in dynamic equilibrium.

4.2.2. Secchi Disk Depth

A Landsat SR median synthetic image for the study year was acquired through the GEE cloud platform; then, the SDD was inverted using the validated model. Finally, the spatial and temporal distributions of the SDD were mapped using the water mask data (Figure 9).

As can be seen from Figure 9, the SDD of the water in Zhejiang was highly dependent on the water depth. In most areas, except for the major rivers and reservoirs, the SDD was less than 70 cm. The downstream SDD of the large rivers was less than 50 cm, while the SDD in the upper reaches was around 100 cm due to the water quality and suspended sediment. On the province-wide scale, the areas with SDDs of 200 cm or greater were the large water bodies located in the mountainous terrain, such as the Xin’anjiang Reservoir in Hangzhou.

To eliminate abnormal SDD values, in this study, the first and last 1% of the data were removed, and then, the data were processed and rounded down to the nearest integer. Then, the water bodies were divided into five categories with equal intervals: low (0–50 cm), medium–low (50–100 cm), medium (100–150 cm), medium–high (150–200 cm), and high (200–267 cm). The data regarding transfers between levels were collected and plotted on a Sankey diagram (Figure 10). The percentage of low SDD (0–50 cm) water bodies was relatively stable; the high SDD water bodies were more obviously variable. Before 1995, there were virtually no areas of water with SDD values of 150 cm or above. During 2000–2010, only a few areas had medium–high SDD (150–200 cm). After 2010, SDD values of greater than 200 cm gradually began to appear.

As can be seen from the Sankey diagram (Figure 7), numerous water bodies with medium–low SDD values transitioned to the medium–high SDD grade after 2010. Between 2010 and 2015, almost 60% of the medium SDD water bodies transitioned to the high and medium–high grades. Additionally, about 50% of the low and medium–low SDD water bodies transitioned to the medium and medium–high grades during this period. Among the different classes of water bodies, those with the largest transformed proportion were the medium–low and medium SDD water bodies, with nearly 50% of these water bodies changing each year, indicating that the stability of the SDD of this class of water bodies was worse.

4.2.3. Suspended Solids Concentration

Using the validated model and a Landsat SR median composite image of the study year, which were obtained based on the GEE cloud platform, the SS concentrations of the water bodies were inverted and a spatial–temporal distribution graph of the concentration was drawn using the water mask data (Figure 11).

As can be seen from Figure 11, the SS concentrations of the water bodies in Zhejiang Province varied considerably between the two poles, and the regional differences were obvious. The SS concentrations in the upstream areas of the large rivers and the various types of reservoirs were less than 20 mg/L. The SS concentrations in the middle and upper reaches of the rivers remained within the range of 40–50 mg/L. The SS concentrations in the lower reaches of the rivers were greater than 200 mg/L, and a small portion of the water bodies within the high-value area had maximum concentrations of greater than 600 mg/L. In the remaining areas, the SS concentrations were mostly in the range of 20–50 mg/L.

Due to the great regional differences in the SS concentration, 40 mg/L was the middle of the largest span of the data, while more than 83% of the data were within 0–60 mg/L (the smallest span of data). Therefore, to distinguish the different concentrations of the water bodies better, in this study, we selected the data with the largest concentration span in the study period, removed the outlier values in the first and last 1%, and classified them according to an interval of 50 mg/L: low (0–10 mg/L), medium–low (10–20 mg/L), medium (20–30 mg/L), medium–high (30–40 mg/L), and high (40–828 mg/L). Then, we investigated the transfers between the different types (Figure 12).

As can be seen from Figure 12, the SS concentrations were concentrated below 20 mg/L and between 40 and 828 mg/L, and the sum of the two concentrations in the middle category accounted for no more than 30% of the total water bodies. Low SS concentration water bodies were virtually absent in 1990 and 1995, and they only had a small distribution from 2000 to 2005. They began to rapidly increase after 2010, accounting for more than 20% of the total water bodies. The percentage of the water bodies with SS concentrations of less than 20 mg/L exhibited a fluctuating yearly increase between 1990 and 2022, with only a slight decrease between 2005 and 2010.

The most significant changes in the SS concentration occurred between 2010 and 2015, during which over 50% of the medium–low SS water bodies transitioned to the low SS category, while nearly 40% of the high SS water bodies transitioned to the medium–low category. The combined proportion of the medium SS and medium–high SS water bodies was less variable and remained in a state of dynamic equilibrium. The reduction in the high SS water bodies was approximately the same as the increase in the combined area of the low and medium–low SS water bodies.

5. Discussion

5.1. Reliability of Research Results

Using permanent water points as a training set, the water areas were identified based on the random forest method and were post-processed using NASA’s SRTM DEM data to obtain water–land classification images. The accuracy validation of the classification results revealed that the kappa coefficients were greater than 0.90, and the overall accuracy (OA) was greater than 0.95% for each study year. Compared with previous surface water extent extraction studies (

Table 5. Comparison of the accuracies of the methods using different satellite data.

Satellite	Research Method	OA (%)	Kappa Coefficient
Landsat 8	RF (this study)	97.36	0.9472
	NDWI	90.81	0.8133
	MNDWI	93.24	0.8610
	AWEIsh	94.52	0.8862
	WI2015	94.48	0.8855

), which employed a variety of different methods such as the NDWI, MNDWI, automatic water body extraction index (AWEIsh), and water identification (WI2015) [17,30,31,32,33], this study applied the EGS operational surface water mapping algorithm, which delivers a better performance in terms of overall accuracy and the kappa coefficient than the four index methods mentioned above [34,35].

In this study of water quality parameter inversion, the inland water body band model utilized in previous studies, in which the Chl-a, SDD, and SS were modeled via regression using Landsat’s Rred, (Rblue + Rgreen), and (Rblue + Rred)/(Rblue/Rred) as the factors, yielded R2 values of 0.846, 0.824, and 0.808, respectively. The inversion results can better reflect the spatial and temporal variations in the water quality parameters in Zhejiang Province and provide a better distinction between regional water bodies.

5.2. Distribution of Surface Water and Water Quality Condition

Hangzhou had a wide area and the largest river and lake in Zhejiang Province, and it possessed the largest surface water area, with an average annual coverage area of 842.65 km². Huzhou and Shaoxing belong to the typical Jiangnan Water Town with a dense network of crisscrossed river courses, and the coverage area was more than 200 km², second only to that in Hangzhou. Quzhou and Lishui in the western part of Zhejiang are located in mountainous and hilly areas. The water system in these areas was underdeveloped, and the area was smaller, only about 120 km². Zhoushan is a coastal archipelago with a narrow geographical area, and the surface water coverage area was less than 10 km².

Table 6. Statistical table of total area and water area change index of each region.

	Area of Districts	Standard Deviation	Area of Water	CV
Hangzhou	16,597	43.78	842.66	0.05
Ningbo	9816	19.72	163.97	0.12
Wenzhou	11,773	7.83	160.42	0.05
Jiaxing	4474	16.53	143.58	0.12
Huzhou	5596	40.06	264.98	0.15
Shaoxing	8778	24.71	227.35	0.11
Jinhua	11,830	21.13	183.68	0.12
Quzhou	8716	28.80	122.72	0.23
Zhoushan	1472	1.10	6.50	0.17
Taizhou	9185	10.61	151.04	0.07
Lishui	17,738	34.82	125.80	0.28

reveals that the CV values of the cities were significantly different. The CV values of Hangzhou and Wenzhou were low (0.05), and the CV values of Lishui and Quzhou were higher (0.28 and 0.23, respectively).

CV is a value with relative meaning. The standard deviation data of Hangzhou is very large, but it has the largest surface water area in Zhejiang Province (large cardinal), which makes its CV index very low. However, the average annual surface water coverage area of Wenzhou is only 160.42 km², and it is mainly concentrated between the two major rivers. Because of the exchange of water bodies between rivers and coastal waters, the water body area presents a relatively stable state. Lishui and Quzhou are located in mountainous and hilly regions. The water system is narrow and long, covering a small area, and there are many small ponds, which are easily affected by climatic factors and cause relatively large area changes. Two large reservoirs were built in Lishui and Quzhou in 2009 and 2001, respectively, with areas of the second and tenth largest in the province (most of the top ten were built before 1990). Therefore, the CV index of these two regions is higher.

The SS concentration was significantly higher in the downstream part of the river in the estuary region than in the other inland water bodies due to the convergence of several rivers and the accumulation of flowing water. Similarly, The Chla was higher and the SDD was lower. Before 2015, the magnitude of the changes in the various types of water in Zhejiang Province was relatively small, but there was a tendency to increase the number of water bodies with better quality conditions. Under the appeal of a policy, the water quality in Zhejiang Province has improved significantly since 2015. In other words, the proportion of the water bodies with high SS and Chla concentrations decreased yearly, the proportion with high Chla concentrations gradually converged to zero, and the proportion with high SSD values gradually increased and reached more than 10%.

5.3. Limitations of the Study

The remote sensing data acquired in this study were synthesized using the median of the Landsat annual satellite data. The data merging methods utilized can affect the original data, thus creating different errors in the study results [36,37,38]. For example, when using the maximum value composite, even if the merged image has been pre-processed by removing clouds, there will still be a large amount of noise in the composite image due to abnormal cloud boundary pixels. When using the minimum value composite, the water region extracted using the random forest method will yield misclassification of mountain shadows, and the most difficult issue to cope with during the post-processing is the misjudgment of the boundaries of the mountain reservoirs and holes in the centers of the water bodies. Therefore, the median value composite of the image was chosen for use in this study, but the median value composite also has the problem of inadequate synthesis of time series image pixels.

The limitation of the ground resolution of the image, the influence of river bridges, and the lack of comprehensive data collected at the sample points resulted in some of the narrow water bodies not being extracted completely, and intermittent patches of water were often formed [39,40,41,42]. In addition, the large scale of the study area, which included a wide variety of water body morphologies in various areas, made it impossible to perform a uniform mathematical morphology operation to fill the entire area.

6. Conclusions

This study explored the extraction of surface water extent based on the EGS operational surface water mapping algorithm and inversion of water quality variables based on validated models using Landsat time series, and created spatial distribution maps of surface water extent and water quality variables in Zhejiang Province at 5 year intervals from 1990 to 2020. The main conclusions of this study are as follows:

(1) The water body extraction accuracy of the EGS operational surface water mapping algorithm, combined with Landsat data, is superior, with kappa coefficients of >0.90 in each of the studied time phases, and the lowest overall accuracies exceed 95.31%. Fragmented mountain water bodies, intricate urban water networks, and meandering rivers can all be extracted well, and the comprehensive mask consisting of multiple indices greatly reduces the potential for the misclassification of land areas.

(2) The inter-annual variations in the surface water bodies in Zhejiang Province exhibited an initially decreasing and then increasing trend, with a minimum value of 2027.49 km² in 2005 and a maximum value of 2610.30 km² in 2015. Based on the statistics on the coefficient of variation (CV) of the water body area in each city in Zhejiang Province, Quzhou and Lishui experienced the most dramatic changes during 1990–2022, and Hangzhou and Wenzhou exhibited the smallest variations in water body area.

(3) According to the multi-spectral data for the water bodies extracted from the results, combined with the water quality band model utilized in previous studies, the water quality parameters in Zhejiang Province were inverted and mapped. In the watershed, the changes in the water quality conditions in Zhejiang Province during 1990–2022 exhibited increases after 2015. The water quality conditions of all of the types of water bodies changed slightly before 2015, while they improved significantly after 2015, and the proportion of the water bodies with better water quality conditions increased.