Remote Estimation of Water Clarity and Suspended Particulate Matter in Qinghai Lake from 2001 to 2020 Using MODIS Images

: Climate change and human activities have been heavily affecting oceanic and inland waters, and it is critical to have a comprehensive understanding of the aquatic optical properties of lakes. Since many key watercolor parameters of Qinghai Lake are not yet available , this paper aims to study the spatial and temporal variations of the water clarity (i.e., Secchi-disk depth, Z SD ) and suspended particulate matter concentration ( C SPM ) in Qinghai Lake from 2001 to 2020 using MODIS images. First, the four atmospheric correction models, including the NIR–SWIR, MUMM, POLYMER, and C2RCC were tested. The NIR–SWIR with decent accuracy in all bands was chosen for the experiment. Then, four existing models for Z SD and six models for C SPM were evaluated. Two semi-analytical models proposed by Lee (2015) and Jiang (2021) were selected for Z SD ( R 2 = 0.74) and C SPM ( R 2 = 0.73), respectively. Finally, the distribution and variation of the Z SD and C SPM were derived over the past 20 years. Overall, the water of Qinghai Lake is quite clear: the monthly mean Z SD is 5.34 ± 1.33 m, and C SPM is 2.05 ± 1.22 mg/L. Further analytical results reveal that the Z SD and C SPM are highly correlated, and the relationship can be formulated with Z SD = 8.072 e − 0.212 C SPM ( R 2 = 0.65). Moreover, turbid water mainly exists along the edge of Qinghai Lake, especially on the northwestern and northeastern shores. The variation in the lakeshore exhibits some irregularity, while the main area of the lake experiences mild water quality deterioration. Statistically, 81.67% of the total area is dominated by constantly increased C SPM , and the area with decreased C SPM occupies 4.56%. There has been distinct seasonal water quality deterioration in the non-frozen period (from May to October). The water quality broadly deteriorated from 2001 to 2008. The year 2008 witnessed a sudden distinct improvement, and after that, the water quality experienced an extremely inconspicuous degradation. This study can ﬁll the gap regarding the long-time monitoring of water clarity and total suspended matter in Qinghai Lake and is expected to provide a scientiﬁc reference for the protection and management of the lake.


Introduction
Lakes can provide fundamental drinkable water resources for human beings and as an important surface water can extensively affect the surrounding climate and ecological environment and vice versa [1][2][3]. With rapid industrialization and urbanization, the water quality of lakes has changed considerably in the past decades [4][5][6]. Increased human activity threatens the water quality of lakes; thus, exploring water quality variation is essential for water security and ecosystem health [7,8]. Qinghai Lake, the largest inland lake and also the largest saline lake in China, is located on the northeastern Tibetan Plateau in Qinghai Province and has a semi-arid and cold climate [9]. Due to its unique geographic location, Qinghai Lake plays a significant role in the climate change of the Tibetan Plateau [10], and climate change on the Tibetan Plateau can heavily affect global climate change [11]. Thus, studies on the ecological environment of Qinghai Lake have drawn much attention in recent years [12][13][14]. Moreover, with the increasing human activities in the Qinghai Lake Basin, the ecosystem of the lake is also being gradually affected [9,15,16]. Hence, exploring how climate change and human activity are affecting the aquatic environment of the lake and how the water quality varies is truly necessary.
Under the dual influence of climate change and human activity, the physical and biochemical properties of Qinghai Lake and its ecological environment have changed significantly. For example, the area of Qinghai Lake shrank before the 2000s and expanded over the past 20 years [12,13,17], and research suggests that precipitation and evaporation were the main driving forces [18]. It was also discovered that the water surface temperature has maintained an upward trend in the last 50 years, and this tendency accelerated in the past few years [12,19]. The increased water surface temperature of the lake was partially affected by global warming. Accordingly, the frozen period of Qinghai Lake was notably decreased as manifested in the delay of the freeze-up and the advance of the ablation [14,20]. With the global concern regarding the carbon cycle of the Earth, the research on carbon sources and sinks of Qinghai Lake was also conducted [21][22][23]. It was reported that the carbon storage of the Qinghai Lake Basin decreased before the 2000s and increased over the past two decades [22], which coincidentally conforms with the water level variation of the lake. Evidence also showed that phytoplankton in Qinghai Lake proliferated in recent decades, boosting the primary production of the ecosystem because of the decreased frozen period and the increased sunshine duration and surface temperature [24]. Generally, recent studies on Qinghai Lake have mainly focused on the variations in water level [17], surface temperature [19], ice phenology [20], carbon dynamics [22], etc. These studies usually tend to use remote sensing images to investigate the spatial distribution and temporal variation of Qinghai Lake on its physical, optical, biochemical, and hydrological characteristics because watercolor remote sensing can provide an economical and convenient way to estimate various properties of the water over a long period.
Remote sensing technology has served as a popular way to estimate water quality parameters, such as chlorophyll (Chl), suspended particulate matter (SPM) and colored dissolved organic matter (CDOM) concentrations, from water-leaving radiance (L w ) or remote sensing reflectance (R rs ) [25,26]. Various inversion models have been developed and applied to inland and coastal lakes from clear to turbid water states [27][28][29][30]. The modeling process usually involves in situ measured water quality properties for model calibration and validation [26]. However, due to the difficult field sampling conditions and the costs of field survey, limited water samples of Qinghai Lake have been acquired and analyzed for a long time-series estimation of the water quality. For example, Shi used the normalized waterleaving radiance derived from the Visible Infrared Imaging Radiometer Suite (VIIRS) to study the chlorophyll-a (Chl-a) concentration and the water diffuse attenuation coefficient at 490 nm (K d (490)) in Qinghai Lake from 2012 to 2018, but it only revealed the variation trend of the lake properties because there is a no actual field measurements of Chl-a and K d (490) for model validation [31]. Generally, the studies on the long-term monitoring of the water quality, such as water clarity (i.e., Secchi-disk depth, Z SD ) and SPM concentration (C SPM ), in Qinghai Lake, are still limited and incomplete.
Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Terra and Aqua satellites can view the entire Earth every 1 to 2 days with spatial resolution from 250 m to 1000 m [32]. This paper took advantage of the MODIS data to derive the Z SD and C SPM quantitatively over the past two decades (from 2001 to 2020) and provide a better understanding of the water quality in Qinghai Lake. To this end, first, the in situ water sampling and field spectral measurement must be collected. Then, the pre-processing routine for remote sensing images needs to be performed to derive the R rs , especially the atmospheric correction (AC). AC is a key step for remote-sensing-based water quality estimation because atmospheric radiation comprises around 90% of satellite-observed radiance at visible wavelengths, and only a small fraction comes from the water body [33,34]. Third, the water quality inversion model should be established with the necessary calibration or validation. At present, empirical, semi-empirical, and semi-analytical models are often employed in remote-sensing-based inversion. The empirical and semi-analytical models are mainly established based on statistical regression between in situ sampled water quality parameters and L w or R rs [35][36][37][38]. Empirical or semi-empirical models are simple but effective given sufficient in situ measurements, but the model transferability is usually limited. The semi-analytical model, in contrast, largely relies on the classic radiative transfer theory and is also accompanied by statistical methods [39][40][41]. Since the semi-analytical model is built based on theoretical physical relations, the estimation accuracy, in practice, may become lower in some complex water bodies with more environmental factors involved.
This paper tested the performance of four commonly-used AC algorithms in Qinghai Lake based on MODIS Level-1 data. The MODIS R rs images were produced with the optimal AC processor. Some empirical, semi-empirical models, and semi-analytical models from works of literature were compared on Z SD and C SPM estimation. Next, the two water parameters in Qinghai Lake were mapped from 2001 to 2020 based on the derived R rs sequences. The relationship between Z SD and C SPM was analyzed using a regression model. In addition, the spatial and temporal variations of the C SPM were studied using time series analysis models. The rest of this article is organized as follows. The study area and data, including field measurements and MODIS images, are introduced in Section 2. Section 3 delineates the details of the AC and time series analysis models employed in this paper. Sections 4 and 5 present the results and discussion, and the conclusion is summarized in Section 6.

Study Area
Qinghai Lake is located in the northeast of Qinghai province in northwestern China. The location of Qinghai Lake (36 • 32 ~37 • 15 N and 99 • 36 ~100 • 47 E) is illustrated in Figure 1. Qinghai Lake is at the highest altitude (~3260 m) of the top largest 50 lakes in the world [31]. The area of the lake is about~4300 km 2 , and it has been getting larger in recent years [12]. The average depth of the lake can reach 21 m [31]. The lake is surrounded by several sub-lakes on the east shore, namely, Haiyan Bay, Shadao Lake, and Gahai Lake, as shown in Figure 1. Qinghai Lake is the largest saline lake in China, and high salinity makes the suspended colloidal particles easily absorb free ions in water to form complexes before settling. Therefore, the lake is often quite clear; moreover, the nutrient level of Qinghai Lake is relatively low [42]. The annual mean precipitation is around 360 mm, and the annual mean temperature is −0.7 • C [24]. It usually enters the freezing period in November and begins thawing in April of the next year [14,20]. In this study, we only investigated the water quality in the ice-free period (i.e., from May to October) in Qinghai Lake.

Field Data
Two field surveys in Qinghai Lake were conducted in September 2008 and August 2016. The sampling points are marked in Figure 1. The sampling points in 2008 were scattered in the center and southern parts of the lake, while the points in 2016 were mainly distributed in the northwestern section. At each station, water samples were collected at depths of around 0 m, 15 m, and 30 m using a Niskin water sampler before quantitative analysis in the laboratory. The spectra were measured using an ASD FieldSpec spectroradiometer following NASA protocols on the above-water radiance measurement [43]. The R rs can be calculated approximately using Equation (1): where S sw , S sky , and S p are the total water-leaving radiance, sky radiance, and the radiance of the gray panel, respectively, measured with the spectroradiometer; the constant r and ρ p are the air-water surface reflectance (a value of 2.5% was used) and the reflectance of the standard reflectance panel, respectively.
There was a total of 29 and 12 valid water samples collected in the surveys of 2008 and 2016. In 2008, only the C SPM was acquired. The SPM was measured by weighing the substance before and after it was filtrated using a Whatman membrane filter with a diameter of 0.45 nm and dried at 105 • C overnight. The SPM was further divided into suspended particulate inorganic matter (SPIM) and suspended particulate organic matter (SPOM) by burning the organic matter out from the SPM at 500 • C for 4 h and reweighing [23]. In 2016, only the water clarity, usually qualified by the Secchi-disk depth, was measured. The Z SD was observed with a 30-cm diameter Secchi disk in the field survey. for the 20-year study. The L1A daily product has a spatial resolution of 250 m, 500 m, and 1000 m for different bands. The L1B product was then generated using the NASA SeaDAS software. Next, four atmospheric correction models were comparatively tested, and the R rs was derived. The generated R rs was automatically resampled to 1000 m for each band. Finally, strict quality control was performed on the R rs product to mask pixels that cannot satisfy the quality standard according to the attached l2_flags metadata of each image. For example, pixels labeled as land, sun glint, cloud, ice, or sensor view zenith angle and solar zenith angle exceeding thresholds were all filtered out. The number of R rs images after processing in each month is distributed as in Figure 2. There is no valid observation in four of the months from 2002 to 2017 due to severe cloud contamination or occasional large-scale missing data in L1A images.

Atmospheric Correction Models
Water is highly absorptive, and the effect of atmospheric gasses and particles must be adequately corrected for further quantitative inversion and analysis [33]. Four atmospheric correction models for ocean color remote sensing, including the NIR and SWIR combined method (NIR-SWIR) [44,45], the management unit of the mathematical models of the North Seas (MUMM) [46], the polynomial-based algorithm applied to the MERIS (POLYMER) [47], and the Case 2 regional coast color processor (C2RCC) [48], were selectively tested in the following experiments. The NIR-SWIR method was established based on the black pixel assumption that the water body shows strong absorption in the NIR and SWIR bands, and the water-leaving radiance is neglectable in that spectral range [49]. For MODIS images, the NIR bands at wavelengths of 748 and 869 nm are used for non-turbid water, and the SWIR bands at 1240 and 2130 nm are for turbid water [45,50]. The MUMM assumes a spatial homogeneity of aerosol properties over a region of interest and replaces the black-pixel assumption with the fact that the ratios of two NIR bands for water-leaving reflectance and aerosol reflectance can be regarded as a constant to build the correction model [46,51]. The additional parameters for the MUMM model, such as the water-leaving reflectance ratio, were all set with the default values in the software. The POLYMER employs a per-pixel spectral optimization method based on a polynomial, and it gains the strength to recover the ocean color in the presence of sun glint [47,52]. The output of POLYMER is the water-leaving reflectance (ρ w ) instead of R rs . The relation between the two can be simply formulated as ρ w (λ) = πR rs (λ) at the wavelength of λ [53]. The C2RCC processor is a set of neural networks trained on simulated top-of-atmosphere reflectance [48,54]. As the name suggests, it is usually applied for Case 2 water. The additional environmental parameters for the C2RCC model were set as follows: the ozone content is 80 DU; the surface air pressure is 801 hPa; the salinity is 12.5 PSU; and the temperature is 20 • C.

Estimation of Z SD
The Secchi disk depth is a common way to measure the clarity of ocean and lake waters. Since limited field data were acquired and the data volume is not sufficient for modeling, we could only test the performance of some existing models and choose the optimal for the following estimation. A semi-analytical model and some empirical or semiempirical models were comparatively tested and verified. Empirical or semi-empirical models proposed by Wu [55], Kratzer [56], and Binding [57] were tested in the following experiments. These models are both implemented based on the R rs at specific wavelengths. For instance, Wu used the blue (469 nm) and red (645 nm) bands to build an exponential model [55]. Kratzer employed similar blue (470 nm) and red (645 nm) bands with a different regression equation [56]. Binding adopted a cubic polynomial regression with one green (550 nm) band [57].
The Z SD derived using the semi-analytical model has been modeled as a function of the downwelling diffuse attenuation and beam attenuation coefficient [58,59]. However, recent studies suggest that the classic theory of water clarity modeling does not really agree with the physical process of the Secchi disk observed by human eyes [39]. The new model can be simply presented as a function of the reciprocal of the diffuse attenuation coefficient at a specific wavelength [39,60]. This paper tested this novel model based on the new theory. The model can be simply formulated as Equation (2): where min(K tr d ) is the K d at the transparent window of the water body, namely the wavelength of maximum transmittance [61] corresponding to the minimum K d for the visible spectral wavelengths between 443 nm and 665 nm. R tr rs denotes the remote sensing reflectance at the wavelength corresponding to the transparent window of the water body. C r denotes the contrast threshold of a human eye when measuring the water clarity with a white Secchi disk, and usually it can be regarded as a constant of 0.013 sr −1 [62]. The K d can be formulated as a function of the total absorption (a) and backscattering coefficients (bb) of water from the perspective of radiative transfer modeling [63]. The parameters of a and bb can be derived using the quasi-analytical algorithm (QAA) [40]. Specifically, the latest version of QAA (QAA_v6) was employed in this paper, and the arithmetic details can be found in [64]. This semi-analytical model has been used in coastal waters [60,65] and in the 50 large lakes on the Yangtze Plain, China [66], and demonstrated its effectiveness in practice [66].

Estimation of C SPM
For the same reason, instead of devising a new algorithm from scratch, some existing models on C SPM were comparatively tested, and the optimal was employed. Typical empirical or semi-empirical models in the literature [53], proposed by Mao [67], Dogliotti [68], Han [69], Novoa [70], and Yu [53], were evaluated in this paper. Among these models, some estimate the C SPM based on the R rs and some on the ρ w with single or multiple spectral bands. Mao used the four customized indices most correlated to C SPM to formulate a proxy variable, and then an exponential equation was established based on this proxy [67]. Dogliotti proposed a generic model based on a single band but differentiated for clear and turbid waters with different feature bands. Moreover, a blending scheme was also designed for waters between the definite clear and turbid states [68]. Han revised Dogliotti's model with different feature bands and an altered blending scheme [69]. Novoa devised a model where the relation between feature bands and C SPM switches among five different condi-tions based on the value of ρ w (671) [70]. Yu established a power-law function based on a proposed generalized index involving various feature bands [53].
The semi-analytical model employed in this paper for C SPM estimation is mainly derived from the particulate backscattering coefficients (b bp ) of water based on the radiative transfer equations [71]. First, the water is classified into four categories: clear, moderately turbid, highly turbid, and extremely turbid, according to the R rs values at the wavelengths of 490, 560, 620, and 754 nm. Second, the radiative transfer equations are used to determine the b bp at the corresponding single wavelength. To be specific, 560 nm for clear water, 665 nm for moderately turbid water, 754 nm for highly turbid water, and 865 nm for extremely turbid water. According to [40], the main process can be formulated as Equation (3): where µ(λ) is a function of R rs at the wavelength of λ, a is the total absorption coefficient, and b bw is the backscattering coefficient of pure water. The calculation of µ and a can follow the procedures specified in the QAA_v6 [64]. The b bw is a known constant corresponding to each spectral band [40,72]. Finally, the C SPM can be estimated by multiplying the b bp with a coefficient (b * bw −1 ) corresponding to each water type at the specific wavelength, as shown in Equation (4): where b * bw denotes the SPM-specific particulate backscattering coefficient. In practice, the b bp -based C SPM model has been widely applied in different water areas with various remote sensing data sources and exhibited its superiority in practice [71,73].

Accuracy Assessment
The model accuracy is assessed and validated by comparing in situ measurements with model estimations in terms of the root mean square error (RMSE), mean absolute percentage error (MAPE), mean normalized bias (MNB), and coefficient of determination (R 2 ). These statistic metrics are calculated as follows: where x i is the i-th in situ measurement, y i is the i-th model prediction,x is the mean of all in situ measurements, and N is the number of valid samples. RMSE and MAPE represent a model's accuracy with regard to prediction errors, and smaller RMSE and MAPE indicate higher accuracy. MNB exhibits a degree of deviation, and the sign of the value suggests the deviation direction. R 2 measures how the predictions can be explained by the model, and a superior model yields R 2 close to one.

Statistical Analysis
The Breaks for Additive Season and Trend (BFAST) and the Mann-Kendall (M-K) test were employed to comprehensively analyze the spatiotemporal variation of Qinghai Lake over 20 years. The BFAST algorithm can decompose time-series data into fitted seasonal, trend, and remainder components for change detection [74,75]. It can be formulated as Equation (9): where Y t , T t , S t , and e t are the observed data value, trend component, seasonal component, and noises at the time t. The model iteratively fits piecewise linear trend and seasonal models to a time series. The magnitude and direction of the changes are derived from the intercept and slope of the components, and then they are used to identify the changes in the time series. Changes occurring in the trend component indicate disturbance, while they demonstrate phenological changes in the seasonal component [75]. The BFAST was used to identify the general trend of water quality variation in Qinghai Lake.
The M-K test can be used for time-series data to ascertain whether there is a consistently increasing or decreasing trend [76]. It is a non-parametric test and works for all distributions. The M-K test is based on a test statistic S under the hypothesis that data are independent and identically distributed [77]. The S can be formulated as Equation (10): where n is the total number of samples, x i and x j are the next and previous samples, respectively, and f is a sign function. A positive value of the normalized S indicates an upward trend, and a negative value suggests a downward trend. Theil-Sen's slope [78] is often used in conjunction with the M-K test to quantify the magnitude of the trend. Considering the seasonal variation in the image series, the seasonal M-K test [79] was employed in this paper to eliminate the effect of seasonality. The M-K test was employed to explore the spatial distribution of the temporal variation trend in a pixel-wise manner. The statistical models were implemented with R programming language in the experiments.

General Properties of Qinghai Lake
Generally, the water of Qinghai Lake was quite clear on the days of sampling. Figure 3a shows the derived R rs curves in the field surveys. Different colors indicate the surveys conducted in different years. Significant reflection peaks at 550 nm, which are usually caused by the weak absorption of chlorophyll, are observed. The R rs in the red and NIR bands (from 700 to 800 nm) are relatively low, which indicates the water was quite clear. Moreover, the in situ measured C SPM and Z SD also confirm that the water quality of Qinghai Lake was in a good condition. The lowest and highest C SPM of the samples are 0.7 and 12.0 mg/L, respectively, and the average is 3.8 mg/L. The SPIM concentration (C SPIM ) accounts for 50.52% of the C SPM regarding all the samples, but the proportion varies considerably, ranging from 23.08% to 85.71%. The remaining 49.48% is from the SPOM concentration (C SPOM ). There are 16 samples in which the C SPM is higher than the C SPIM among the 29 samples. The Pearson correlation coefficient (r) between the C SPM and the C SPIM was as high as 0.91, and the relation approximates an exponential function with Equation (11) (RMSE = 0.66, MAPE = 34.78%, MNB = −7.65%, R 2 = 0.81): The fit curve is shown in Figure 3b. The symbol N indicates the valid samples for evaluation. The red line is the regression line, and the bandwidth in translucent green indicates a 95% confidence interval. With respect to the in situ measured Z SD , the shallowest and deepest of the Z SD are 1.4 and 5.1 m, respectively, and the average is 3.13 m among the 12 valid samples.

Accuracy of Atmospheric Correction
The performance of atmospheric correction models was evaluated by comparing the in situ measured R rs with the remote sensing derived R rs in a point-to-point manner.
The point values extracted from remote sensing images corresponding to in situ stations were averaged with a 3 × 3 neighboring window to reduce accidental errors. Because of the severe cloud contamination in the remote sensing images on some days of the filed surveys, only 14 valid matched points were obtained after data filtering. The scatter plots in Figure 4 exhibit the fitness of the four atmospheric correction models on each spectral band. The quantitative metrics of each model are presented in Table 1. The best ones are marked in bold. There is no result produced at the wavelengths of 469 nm, 555 nm, and 645 nm with C2RCC and POLYMER; hence, the corresponding cells remain blank. Several conclusions can be drawn from Figure 4 and Table 1. First, the C2RCC produces the highest uncertainty in all spectral bands (MAPE ≥ 66%, MNB ≤ −210%). The C2RCC corrected R rs is significantly lower than the in situ measurement. Second, the POLYMER functioned well at the blue and green wavelengths (MAPE ≤ 14%, |MNB| ≤ 6%, R 2 ≥ 0.64), but it showed poor performance at the red wavelengths (MAPE ≥ 24%, |MNB| ≥ 33%, R 2 ≤ 0.50). Third, the NIR-SWIR and MUMM were nearly comparable. The NUMM obtains higher accuracy from the blue to red bands (RMSE ≤ 0.003 sr −1 , MAPE ≤ 30%, |MNB| ≤ 6%, R 2 ≥ 0.45), but the MUMM performed defectively at the wavelength of 748 nm with a significantly enlarged outcome (MAPE = 183.21%, MNB = 58.42%, R 2 = 0.07). The accuracy of NIR-SWIR (RMSE ≤ 0.004 sr −1 , MAPE ≤ 41%, |MNB| ≤ 92%, R 2 ≥ 0.34) is slightly lower than MUMM, but it produces averagely acceptable results for all spectral bands. Fourth, the results with NIR-SWIR (MNB less than 0 for 10 bands out of 11) and MUMM (MNB greater than 0 for 9 bands out of 11) were generally a little lower and higher than the in situ measurements, respectively.
Band ratios are often used to measure the system errors induced by the atmospheric correction [5], and some band ratios are also illustrated in Figure 4. The ratios of 443/547 and 488/547 were close to the identity line with the NIR-SWIR, MUMM, and POLYMER. The ratios of 667/547 with the NIR-SWIR and MUMM were slightly lower and higher than the identity line. The NIR-SWIR achieved the best at the ratio of 748/547. In conclusion, the NIR-SWIR and MUMM performed better than the other two models, but considering the MUMM did not perform well in the NIR band (748 nm), and the NIR band is often essential in estimation models of water-quality parameters, the NIR-SWIR was chosen to conduct the atmospheric correction in the following experiment.

Estimation of Z SD and C SPM
Several empirical, semi-empirical, and semi-analytical models were tested for the Z SD and C SPM estimation. The model-estimated results were compared with the in situ measurements for validation. The quantitative evaluation of the four models for Z SD and the six models for C SPM are listed in Tables 2 and 3, respectively, in terms of RMSE, MAPE, MNB, and R 2 . Generally speaking, the empirical and semi-empirical models showed a dreadfully poor performance concerning all the statistical metrics (RMSE ≥ 0.7 m, MAPE ≥ 26%, |MNB| ≥ 36%, R 2 ≤ 0.08 for Z SD ; RMSE ≥ 2.6 mg/L, MAPE ≥ 46%, |MNB| ≥ 21%,  Figure 5. The red line is the regression line, and the translucent green area is the 95% confidence interval. The gray line is the 1:1 line for reference. The R 2 can reach 0.7, and the RMSE is less than 0.8 for both the Z SD and C SPM with the limited samples. Basically, the MOSIS-derived Z SD and C SPM are slightly lower and higher than the in situ measurements, respectively. The points in Figure 5a for Z SD are distributed comparatively centralized, while they are more dispersed in Figure 5b for C SPM with a wilder confidence interval. Since there are only 11 valid samples for C SPM , we believe the results are acceptable. Hence, the semi-analytical models were adopted for the following long-series water quality estimation.   Figure 5. Comparison between the in situ measurements and MODIS-derived Z SD (a) and C SPM (b) with semi-analytical models (the red line is the regression line, and the translucent green indicates a 95% confidence interval).
The daily Z SD and C SPM in non-frozen periods (from May to October) were produced with the MODIS R rs data from 2001 to 2020. Figure 6 illustrates the mean variations of Z SD and C SPM calculated by averaging pixel values of the whole lake on the daily and monthly scales. Statistically, the monthly mean Z SD is 5.34 ± 1.33 m, and C SPM is 2.05 ± 1.22 mg/L. The variation of both Z SD and C SPM exhibits an observable seasonal trend in Figure 6. The Z SD continues decreasing and C SPM keeps increasing throughout the year. Z SD and C SPM are highly correlated parameters manifested in the fact that higher C SPM usually leads to lower Z SD and vice versa. A regression analysis was performed on the monthly averaged Z SD and C SPM to identify the relationship between their quantality. Figure 7 reveals the two are exponentially correlated with Equation (12) (RMSE = 0.25, MAPE = 4.01%, MNB = −0.8%, R 2 = 0.65): The density plot in Figure 7 illustrates the distribution of the corresponding pixel values in the Z SD and C SPM images. The color indicates the number of pixels in the hexagon bin. The red line is the regression line with a satisfactory fitness. Since the relationship between Z SD and C SPM has been qualified, only the C SPM is discussed regarding its spatial and temporal variation for the sake of brevity.

Spatial and Temporal Variation of C SPM
Distinct seasonal cycles can be observed in Figure 6. To further explore the seasonal variation, Figure 8 exhibits the monthly spatial variation of C SPM . Each of the subfigures demonstrates the spatial distribution of the C SPM where the pixel values were averaged for every corresponding month over the 20 years. The C SPM continues to increase from May to October in a year. The turbid water was mainly distributed along the lake edge every month, specifically on the northwestern and southeastern shores. Because Haiyan Bay is so close to the main Qinghai Lake, there are distinct mixed pixels, sometimes referred to as the land adjacency effect [80], leading to poor water quality on the southeastern shore in Figure 8.   , the water quality varied somewhat, and there was not much change observed during this period. In 2020, the lake again experienced a trivial amelioration in water quality. Moreover, the BFAST algorithm was employed on the monthly mean C SPM to explore the seasonal and total trends quantitatively. The result is revealed in Figure 10. The total variation in the month scale is fairly complicated. By decomposing it, the total trend over the 20 years is that the water turned less clear distinctly from 2001 to 2007; it gained sudden distinct recovery in 2008, and after that, the C SPM in the lake has been gradually increasing. By removing noises, the total trend is basically similar to the above conclusion observed in Figure 9. Significant seasonal degradation is observed in Figure 10, confirming the conclusion from Figures 6 and 8. However, the variation of segregated noises is substantially volatile, which shows that temporal change is a complicated process despite the overall decreasing trend.
A seasonal M-K test was employed to determine the spatial variation trend from 2001 to 2020 with the monthly mean data of every year. Figure 11 presents the spatial distribution of the M-K test result and the corresponding Theil-Sen's Slope to identify the magnitude of the C SPM variation trend. The hypothesis that there is a continuous increasing or decreasing trend is accepted with a significance level of 5%, and the trend is determined to be significant with a significance level of 1%. Generally, the main body of the lake tended to become less clear regarding C SPM , and the area showing an increasing trend accounted for 81.67% of the total. The water body near the shore where the water is usually more turbid than the main body showed no trend or had some improvement. The area with decreasing turbidity only accounted for 4.56%. The slope in Figure 11b is generally consistent with the M-K trend in Figure 11a. The biggest magnitude of the trend occurred on the northwestern and southeastern shores regarding the Theil-Sen's Slope in Figure 11b.

Discussion
AC is a necessary preprocessing step for watercolor remote sensing applications due to the high absorption of water. Since the water in Qinghai Lake is less turbid, some AC algorithms, including NIR-SWIR, MUMM, C2RCC, and POLYMER, designed for oceanic applications over open and coastal waters were tested in this paper. The results of NIR-SWIR were constantly lower than the in situ measurements in the red bands. This might be because the black pixel assumption may not always be valid [51], and the extrapolation process on atmospheric influence from NIR to visible bands might cause some errors to lower the estimation. The MUMM performs the best on nearly all the bands, but it failed for the band at the wavelength of 748 nm with a much higher estimation. This might be because, first, the signals of NIR spectra are extremely weak in Qinghai Lake; second, the in situ average R rs (748)/R rs (869) is 1.870 ± 1.740 and the deviation is rather large, not fully satisfying the underlying assumption of MUMM. The POLYMER can compete with the MUMM in the blue and green bands but fails in the red and NIR bands. This malfunction might be ascribable to the low reflectance in the NIR bands. The C2RCC produced the highest uncertainty in Qinghai Lake. Considering the C2RCC is a machine-learningbased model, its transferability in different regions cannot be guaranteed, and further model tuning may be needed. Moreover, the C2RCC needs additional environmental factors as inputs, and the estimated values over the study period are set empirically in the experiments. These less accurate parameters may lead to some unexpected errors.
After assessing various AC algorithms, some available water quality inversion models were evaluated. The water clarity and C SPM are important parameters for water quality monitoring, so the two water parameters were evaluated using inversion models with the acquired in situ water samples. Generally, the semi-analytical models achieve the best performance for the Z SD and C SPM estimation, while the existing empirical and semi-empirical models are not suitable in Qinghai Lake. Nevertheless, this is understandable because empirical and semi-empirical models were established using water samples gathered in a specific region. An empirical model may only fit the relation well in that specific region and does not possess favorable transferability. The semi-analytical model, in contrast, can successfully handle the situation in Qinghai Lake where the water is less turbid. The estimation results for Z SD and C SPM are slightly lower and higher than the in situ measurements. This may be because there are external errors introduced from the AC process, and the semi-empirical model is completely uncalibrated in Qinghai Lake due to the lack of sufficient in situ data. However, overall, the model accuracy is satisfactory for both Z SD and C SPM with R 2 > 0.7. The Z SD and C SPM are highly correlated because the phytoplankton abundance in Qinghai Lake is not high [81], and the water clarity is mainly affected by the C SPM . The water clarity and C SPM estimated in this paper can only manifest partial optical and chemical characteristics of the water. Some biochemical properties, such as Chl-a and CDOM, still need further investigation. The contribution of this paper is that it demonstrated a feasible way to conduct long-term water quality monitoring with MODIS images in Qinghai Lake, and it can be applied to other watercolor parameters.
By further analyzing the derived C SPM series, the seasonal deterioration of water quality can be found. A similar trend was also disclosed regarding the Chl-a and K d (490) in Qinghai Lake [31]. This may be because the lake is often frozen from November to the April, and it is hardly affected by the external environment. Meanwhile, the self-purification of the lake can make the water even more clear. After the thaw in April, the lake is gradually altered by climate changes and human activities. According to studies on the topic [82], we speculate that this may be involved with land use and climate change, but it still remains open to explore. The annual variation is complicated, but the overall trend is that the C SPM has been rising over the 20 years, except for the sudden distinct decline in 2008. The government of Qinghai province implemented a plan in 2008 aimed at ecological environment protection and comprehensive management of the Qinghai Lake Basin [83], which may contribute to the water quality improvement in that year. The water quality was degraded significantly before 2008, and it was still experiencing degradation after 2008, but the change was quite slow. The spatial variation mainly happened along the lakeshore, especially in the northwest and the southeast. There were no significant upward or downward trends in the edge areas. Generally, the main body of Qinghai Lake witnessed increased C SPM and decreased Z SD over the 20 years, but the changes were at an extremely slow and smooth pace. While few studies have looked into the long-term monitoring of Qinghai Lake over recent years, the work in this paper can provide a preliminary estimation of the water quality variation.

Conclusions
This paper devotes itself to exploring the spatiotemporal distribution and variation of Qinghai Lake from 2001 to 2020 using MODIS images. First, four atmospheric correction models including NIR-SWIR, MUMM, POLYMER, and C2RCC were evaluated. Generally, the NIR-SWIR and MUMM outperform the POLYMER and C2RCC; the performance of NIR-SWIR and MUMM is comparable: the MUMM shows higher accuracy in blue and green bands, but the NIR-SWIR gains higher accuracy in the NIR bands. Then, the NIR-SWIR algorithm was employed for the atmospheric correction of the MODIS time series. Second, several existing models for Z SD and C SPM estimation were tested, and the results revealed that the semi-analytical models are superior to the empirical models in Qinghai Lake. Third, he semi-analytical models proposed by  and Jiang (2021) were adopted to map the distribution of daily Z SD and C SPM in Qinghai Lake over the 20 years. Fourth, the quantitative relationship between Z SD and C SPM was established, and the BFAST analysis and M-K test were performed to understand the spatiotemporal variation of the C SPM during the years.
Several conclusions can be drawn from the distribution and variation of the water quality. First, the C SPM and C SPIM and the Z SD and C SPM are highly correlated water quality parameters: the C SPM shows an exponential relation with C SPIM , and the Z SD is an exponential function of C SPM . Second, turbid water mainly existed along the edge of Qinghai Lake, especially on the northwestern and northeastern shores. Third, substantial water quality deterioration can be observed seasonally. In other words, the water became constantly less clear from May to October annually. Fourth, the water quality was declining generally, but the variation was complex. Specifically, the water quality broadly deteriorated from 2001 to 2008; 2008 witnessed a sudden distinct improvement, and after that, the water quality experienced an inconspicuous degradation. Fifth, there is 81.67% of the total area with constantly increased C SPM , and the area with decreased C SPM occupies 4.56% of the total. This study can fill the gap regarding the long-term monitoring of water clarity and total suspended matter in Qinghai Lake, and hopefully, it can provide a scientific reference for the protection and management of Qinghai Lake.