Geographic Characteristics and Meteorological Factors Dominate the Variation of Chlorophyll‐a in Lakes and Reservoirs With Higher TP Concentrations

Nutrients such as phosphorus and nitrogen lead to extensive growth of harmful algae in lakes and reservoirs, which results in eutrophication. The driving mechanism of primary productivity change in lakes and reservoirs at a wide spatial and temporal scale remains largely unknown. We establish a water quality database using a stacking machine learning model, including monthly chlorophyll‐a (Chl‐a), total nitrogen (TN) and phosphorus (TP) concentrations in 255 lakes and 332 reservoirs from 1980 to 2018. We find an increase in the number of lakes and reservoirs at risk of algal blooms, with approximately 2.66% exhibiting Chl‐a concentrations exceeding 30 μg L−1 by 2018. Variations in Chl‐a concentrations were not entirely synchronized at the individual and regional levels with TN and TP concentrations, or their stoichiometric ratios, as estimated by a hierarchical linear model spanning 1980 to 2018. Further, we discovered unexpected effects of geographical features (i.e., latitude, longitude, and slope) and meteorological factors (i.e., air temperature) on Chl‐a concentrations in the lakes and reservoirs with higher TP concentrations (>0.041 mg L−1 for lakes and >0.027 mg L−1 for reservoirs), through the use of multiple regression trees and structural equation model analysis. Our findings underscore the importance of (a) implementing flexible nutrient pollution control approaches based on nutrient ecoregions that consider geographical variations, and (b) developing mitigation and adaptation strategies to address the uncertain risks posed by climate change in the prevention and control of eutrophication in lakes and reservoirs.


Introduction
Primary productivity, the fundamental function of aquatic ecosystems, is influenced by the variation of phytoplankton biomass, making it a crucial topic in limnology (Liang et al., 2020;Oliver et al., 2017;Olson & Jones, 2022;Paltsev & Creed, 2022;Peng et al., 2021;Wu et al., 2022;Xu et al., 2021).Global satellite remote sensing data indicate that phytoplankton blooms in large lakes have increased since the 1980s (Ho et al., 2019;Hou et al., 2022).The production of algal toxins, odorous substances, surface scum, and decreased dissolved oxygen by these freshwater blooms adversely affect aquatic organisms, recreation and tourism, drinking water supplies, and irrigation (Huisman et al., 2018;O'Neil et al., 2012).In the United States alone, freshwater blooms result in annual economic losses of $4 billion, and China is no exception (Huisman et al., 2018;Xu et al., 2021).Resource constraint theory has been employed to identify the limiting factors for the growth in phytoplankton biomass, with nitrogen (N) and phosphorus (P) widely recognized as key limitations (Conley et al., 2009;Merz et al., 2023;Paerl et al., 2016).
Remedial measures have been initiated in developed countries since the mid-twentieth century to combat lake eutrophication, including source control and pollution interception (Ho et al., 2019;Hou et al., 2022;Huisman et al., 2018).China's rapid development following the reform and opening-up policies of 1980 led to the implementation of total pollutant control around 2000 and intensified efforts to prevent eutrophication after 2012 (Huang et al., 2021;Ma et al., 2020).Under the balance of economic development and environmental protection in China in the past 40 years, understanding the spatio-temporal variation trends and driving mechanism of phytoplankton biomass in China's lakes and reservoirs can provide a framework for studying the relationship between nutrients and phytoplankton, offering insights for lake and reservoir protection and eutrophication control in developing countries (Wang et al., 2022).
While several studies have explored the empirical relationships between the concentrations of nutrients and chlorophyll-a (Chl-a) and the underlying mechanisms influencing phytoplankton biomass, our knowledge of these factors remains incomplete (Davidson et al., 2023;Olson & Jones, 2022;Paltsev & Creed, 2022;Peng et al., 2021;Quinlan et al., 2020), in part, because of three predominant barriers.First, there is a lack of long-term and large-scale systematical data that can be used to reveal the temporal and spatial patterns in phytoplankton biomass (Ma et al., 2023).Historically, two methods have primarily been used to characterize phytoplankton biomass through time in lakes and reservoirs, including analyzing remotely sensed images and measured surface water Chl-a.While remote sensing offers wide coverage, low cost, and long-term monitoring, it is plagued by unfavorable observational weather conditions at the time of collection (Ho et al., 2019).In contrast, direct measurement of Chl-a provides more accurate information on phytoplankton biomass (Huo et al., 2019), but is more difficult to consistently collect over broad areas.China initiated the construction of a water environmental monitoring network in 1989, but large-scale standardized monitoring data were lacking during its early stages (Huang et al., 2021).Second, studies on the temporal and spatial changes of phytoplankton biomass and/or nutrient status in large-scale lakes have often overlooked multi-level nested relationships, such as nested timelacustrine-basin relationships (Oliver et al., 2017).Temporal and spatial variations in Chl-a concentration are commonly determined through linear or nonlinear fitting at the reginal scale, which disregards the random influences at lake scales (Filstrup & Downing, 2017).Third, investigations into the driving mechanisms of phytoplankton biomass have primarily focused on the relationships between nutrients and algae.Empirical models have revealed that the relationship between TP and Chl-a typically follows a sigmoidal function, and some studies have noted variations in the TP-Chl-a relationship with latitude (Olson & Jones, 2022;Quinlan et al., 2020).However, it is crucial to consider the specific ecosystem when examining the relationship between nutrients and Chl-a.Research has shown that P control is effective for low eutrophication or deep-water lakes, but may not be suitable for shallow water lakes with high eutrophication due to endogenous cycling (Paerl et al., 2016;Qin et al., 2020).
To address these knowledge gaps, this study aims to: (a) simulate monthly Chl-a concentrations in 255 lakes and 332 reservoirs in China from 1980 to 2018 using a stacking machine learning model, (b) analyze the spatiotemporal variations of Chl-a, TN, TP, and the TN:TP mass ratio in these lakes and reservoirs between 1980 and 2018 through a hierarchical linear model (HLM) that considers time fixed effects as well as lake-(or reservoir-) and region-specific random effects, and (c) determine the driving mechanism behind the variations in Chl-a concentration using multiple regression tree (MRT) and structural equation model (SEM) analyses.This study attempts to provide the "big-picture" view of the long-term historical trends in trophic status and primary productivity in China's lakes and reservoirs.

Investigative Framework
The general architecture of the study consists of three main steps, namely simulation, variation, and attribution of Chl-a concentrations (Figure 1).In the first step, the relationships between Chl-a concentration and the predictor variables including aquatic TN and TP concentration, along with meteorological, topographical, soil, and landuse characteristics, as well as social and economic status, were established and predicted using a two-layer learning approach.We used a stacked machine learning model to predict 389,181 Chl-a concentrations (including 274,716 monthly, 91,572 seasonal and 22,893 annual values) in 255 lakes and 332 reservoirs from 1980 to 2018.The model was based on 39,466 Chl-a measurements at monitoring sites between 1990 and 2018.In the second step, a HLM including time fixed effects and lake-and region-specific random effects was used to estimate temporal trends in Chl-a concentrations from 1980 to 2018.Lakes and reservoirs with increase, decrease, or no change in Chl-a concentration were classified based on whether the associated 90% confidence interval of the lake-or reservoir-specific trend in Chl-a concentration overlapped zero.In the third step, MRT and SEM analyses were performed to identify landscape, climate and nutrient controls associated with Chl-a concentrations.

Data Sources and Preprocessing
Observed water quality data of Chl-a concentration for 255 lakes and 332 reservoirs from 1990 to 2018 were provided by the China National Environmental Monitoring Center, Provincial Environmental Monitoring Center and Research Institutes (RI) (Figure S1 in Supporting Information S1).The concentrations were determined using the spectrophotometric method recommended by the China National Environmental Protection Agency (CNEPA, 2002).Monthly aquatic TN and TP concentrations in all 255 lakes and 332 reservoirs were obtained from our previous study (Ma et al., 2023).The geographic location for each lake or reservoir was downloaded from the China Lake Scientific Database (http://www.lakesci.csdb.cn/).The meteorological data were collected from the CN05.1 data set produced by the China Meteorological Administration with a horizontal resolution of 0.25°× 0.25° (Wu & Gao, 2013).The land-use data set of China from the period of 1980-2018 at a resolution of 30 ✕ 30 m was acquired from the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences.The average net anthropogenic N/P input (NANI/NAPI) for each watershed was derived from the relationship between population density (POP) and NANI/NAPI (Ma et al., 2023).Spatial data of Gross Domestic Product and POP at a resolution if 1 × 1 km were obtained from the Resource and Environment Science and Data Center (https://www.resdc.cn/)(Liu et al., 2005).Soil properties were extracted from the 1:1,000,000 scale digital soil properties map of the Institute of Soil Science, Chinese Academy of Sciences (http://vdb3.soil.csdb.cn/).The observed Chl-a data and the predictor variables were transformed, and normalized to eliminate the effects of measurement differences that range in the orders of magnitude.A total of 54 candidate variables (Table S1 in Supporting Information S1), including nutrients, geomorphology, climate, land use, soil properties and socioeconomic data were used for variable filtering and model building for Chl-a simulation.For more details on data sources, transformation, and cleaning, see the Text S1 and S2 in Supporting Information S1.

Machine Learning Algorithms
Input variable selection.The important input variables influencing Chl-a concentrations were identified by the maximal information coefficient (MIC) method (Reshef et al., 2011).After min-max normalization of the data, the correlated relationships between candidate variables and log Chl-a concentration were estimated using MIC scores.Variables with MIC scores larger than 0.2 were selected as predictor variables for the Chl-a simulation (Reshef et al., 2011).Furthermore, to avoid multicollinearity and ensure the independence of predictor variables, variables with Spearman correlation coefficients larger than 0.9 were excluded, whereas variables with significant correlated relationships with Chl-a concentration were retained (Ma et al., 2023).
Base models.Due to the nonlinear, non-normal, and heterogeneous relationships between Chl-a concentration and predictor variables, three nonlinear based models were used, including K-Nearest Neighbors (KNN), Random Forest (RF) and Support Vector Machine (SVM).These models were selected due to their high popularity and the performance achieved during previous studies (Almuhtaram et al., 2021;Olson & Cormier, 2019;Yan et al., 2019).KNN, a lazy learning algorithm, uses the average value of the sample output of the nearest k samples as the regression value (Zhang & Zhou, 2007).The number of neighbors considered, k, was optimized to 5 by evaluating the coefficient of determination (R2 ).RF regression constructs randomized and independent decision trees for each bootstrap sampling subset and obtains simulation results by parallel methods (Olson & Cormier, 2019).Two parameters in the RF model, the number of trees to grow (ntrees) and the number of variables randomly sampled for splitting at each node (mtries), were optimized to 500 and 2, respectively, using a grid search hyperparameter optimization.The SVM regression algorithm is mainly based on insensitive and kernel function algorithms (Zheng et al., 2011).When the fitted mathematical model is used to express a curve in the multidimensional space, the result obtained by the insensitive function is the "pipeline" containing the curve and the training point.Of all the sample points, only those on the "wall" determine the location of the pipe, which is called the "support vector."To adapt to the non-linearity of the training sample set, a kernel function was introduced for dimensional enhancement and to control overfitting.The regularization parameter (C) and the hyperparameter (sigma) of the radial basis kernel functions were optimized to 2 and 1, respectively, for the best performance of the SVM model.The KNN, RF, and SVM models were built using the "class," "randomForest," and "kernlab" packages, respectively, in the R software package.Default parameters were used with the exception of the aforementioned parameters (Karatzoglou et al., 2022;Liaw & Wiener, 2002;Venables & Ripley, 2002).
Model stacking.A two-layer stacking modeling was used here, in which a higher-level model integrated the results of lower-level models to improve the predictive power of the regression (Ghahramani, 2015).The generalized linear model was used to combine the outputs of the three base models (KNN, RF, and SVM) (Wang et al., 2021).The simulation function of the linear stacking model was as follows: where y is the stacking target (i.e., the final simulation of Chl-a concentrations), and w i and f i represent the base model weight and prediction from M individual algorithms (M = 3 in this study).The optimal set of weights was estimated by quadratic programming, as follows: where y i is the observed value of the ith observation, x i is the ith observation composed of all environmental variables, f i m ( x i ) is the simulation of the ith simulation during the base model m training with the ith observation removed data set, and ω st is the objective function which is estimated by the least square linear regression of y i on the linear combination of f i m ( x i ) under two constraints.The R package routine "caretEnsemble" was used to assign the stacking model (https://github.com/zachmayer/caretEnsemble).
Model validation and evaluation.The 10-fold cross-validation (CV) method was used to optimize the hyperparameters and to evaluate the predictive and hindcast performance of the model.The observation data set was randomly split into 10 subsets, each containing 10% of the data.In each turn, 9 subsets were used to train and the remaining one was used for model validation.The performance of the CV simulations was estimated by four statistical metrics, including R 2 , Nash-Sutcliffe efficiency (NSE), root mean square error (RMSE), and mean absolute error of the linear regression between the observed and predicted Chl-a concentrations (Legates & McCabe, 1999).More details of model evaluation are provided in the Text S3 in Supporting Information S1.

Hierarchical Linear Models
Based on the above predicted data set of monthly Chl-a concentrations for 255 lakes and 332 reservoirs from 1980 to 2018, nutrient trends were examined using HLMs.The heterogeneity of local and regional drivers was expected to form a natural sequence of nested and hierarchical structures in the data, where observations within a lake might be more similar than those across lakes, and lakes within regions might be more similar than lakes across regions (Oliver et al., 2017).The variation within and across lakes and regions was estimated in the HLMs by examining the relationship between sample year and response (i.e., the "trend"), which differed by lake and region (i.e., lakeand region-specific trends or "random effects").Annual TN, TP, TN/TP ratio, and Chl-a concentration data from the lakes and reservoirs were used for the following HLM analysis.In doing so, the first level of hydrological units for China was used at the region level of the model.See more details of the model in the Text S4 in Supporting Information S1.

Statistical Analysis
MRT and SEM analyses were used to determine which geographic, climatic, and nutrient variables explain the variation in Chl-a concentration in lakes and reservoirs.MRT takes independent variables as classification nodes and uses the binary split method to progressively divide dependent variables into categories that are as homogeneous as possible.Cross-validation and the cost complexity parameter were used to prune the branch to obtain the regression tree with minimized error and appropriate tree size (De'Ath, 2002).Compared to MRT, SEM can identify the indirect effects between two explanatory variables.An a priori SEM was applied to assess the causal relationship between the predictors and Chl-a concentration, describing the expected relationship between variables.The probability level of the chi-square test (P > 0.05), the root mean square residual (0 ≤ RMR ≤ 0.05), the comparative fit index (0.90 ≤ CFI ≤ 1.00), and the RMSE of approximation (0 ≤ RMSEA ≤ 0.05) were used to assess the goodness of fit and determine the best performing model (Arhonditsis et al., 2006).The R packages "rpart," "lavaan," and "semPlot" were used to perform the models (https://cran.r-project.org/).

Temporal Variation in Nutrients and Chl-a in Lakes and Reservoirs Since 1980
The simulation of Chl-a concentration has high accuracy using a stacking model with an ensemble data set including annual, quarterly, and monthly data (see more details in the Text S5 in Supporting Information S1).The predicted data set revealed an increasing trend in Chl-a, indicating a growing threat of algal blooms to lakes and reservoirs in China from 1980 to 2018 (Figure 2).The proportion of lakes and reservoirs with annual mean Chl-a concentrations below 10 μg L 1 decreased from 68.50% to 89.46% in 1980 to 51.97% and 79.82% in 2018, respectively.Moreover, the proportion of lakes and reservoirs with annual mean Chl-a concentrations greater than 30 μg L 1 increased from 0.39% to 0% in 1980 to 2.36% and 0.30% in 2018, respectively.To determine the causal In general, lake Chl-a trends were positively correlated with changes in TN and TP, whereas Chl-a in reservoirs responded primarily to changes in TP (Figures 4g and 4h).When considering the confidence intervals for rates of change, significant changes in TN and Chl-a were observed in 19% of lakes and reservoirs, while no significant changes were found in 41% of lakes and 44% of reservoirs (Figures 4c and 4d).Regarding TP and Chl-a, 35% of lakes and 33% of reservoirs showed no significant changes, while 23% of lakes and 26% of reservoirs exhibited  Water Resources Research 10.1029/2023WR036587 significant changes in both TP and Chl-a (Figures 4e and 4f).Notably, points in the fourth quadrant of Figure 4 indicate significant changes in both TP and Chl-a, suggesting that an increase in TP concentrations in lakes or reservoirs is unlikely to lead to a decrease in Chl-a (Figures 4e and 4f).The contour plot also reveals vertical bands, indicating that Chl-a primarily responds to changes in TP (Figures 4g and 4h).Interestingly, in reservoirs, the red area in the third quadrant of the contour plot indicates nine reservoirs where TN and TP concentrations decreased, but Chl-a increased, such as within the Dongpinghu, Liujiaxia, and Qinshan reservoirs (Figure 4h).In summary, the lack of complete synchronization in TN, TP, and Chl-a concentrations from 1980 to 2018 at both individual and regional levels of lakes and reservoirs in China suggests that the relationship between nutrients and Chl-a, inherent in the commonly used traditional paradigm, may not be universal (Davidson et al., 2023;Oliver et al., 2017;Quinlan et al., 2020).These findings emphasize the complex and context-dependent nature of nutrient-Chl-a dynamics in aquatic ecosystems.
Changes in the concentrations of TN, TP and Chl-a in lakes and reservoirs were not completely synchronous between 1980 and 2018 (Figure 3, Table 1).On average, TP in lakes and reservoirs exhibited a negative trend, with a decrease of 0.10% year 1 and 0.52% year 1 from 1980 to 2018, respectively, while the stoichiometric ratio (TN:TP) showed a positive trend, increasing at a rate of 0.26% year 1 and 0.30% year 1 , respectively (Table 1).However, at the individual lake and reservoir scale, the trends in TP concentration and TN:TP were more diverse.Approximately 28% of lakes and 43% of reservoirs displayed a decreasing trend in TP concentration, while 21% of lakes and 6% of reservoirs showed an increasing trend.In contrast, only 12% of lakes and 8% of reservoirs exhibited a decreasing trend in TN:TP, while 29% of lakes and 23% of reservoirs showed an increasing trend.The imbalance in stoichiometric variation between lakes and reservoirs was evident in the relationship between the change in TN and TP (Figures 4a and 4b) and may be related to the rapid improvement in municipal wastewater treatment (Tong et al., 2020).
The changes in Chl-a concentration were not fully synchronized with the variations in TN, TP, and TN:TP at the regional level (Figure 3, Table 1).For instance, TN and TN:TP of reservoirs in the Haihe River Basin decreased, while Chl-a increased.However, in the Huaihe River Basin, TP, TN:TP, and Chl-a in lakes all decreased over time.Only about 20% of the regions exhibited changes in TN, TP, or Chl-a within the lakes or reservoirs.Additionally, the distribution of lakes with elevated Chl-a concentrations was primarily observed in the middle  and lower reaches of the Yangtze River, the Yunnan-Guizhou Plateau, and the lower reaches of the Yellow River basin.Regions with increases in Chl-a were less contiguous (Figure 3).Moreover, the population means of TN, TP, Chl-a, and TN:TP concentrations from 1980 to 2018 were generally not correlated with the estimated rates of change in lakes or reservoirs (Figure S2 in Supporting Information S1).These findings highlight the complex dynamics and spatial heterogeneity in the relationships between nutrient concentrations and Chl-a variations in lakes and reservoirs at both the population and regional levels.

Influence of Geographical Characteristics and Meteorological Factors on Chl-a
To further investigate the variations in Chl-a concentrations and their attributions, we conducted an analysis using MRT and SEM (Figures 5 and 6).Our findings indicate that TP concentration played a crucial role in predicting Chl-a concentration, acting as the primary predictor in the regression tree.The regression tree's first split occurred at 0.041 mg L 1 in lakes and 0.027 mg L 1 in reservoirs, clearly distinguishing the two ecosystem types.Lakes exhibited higher TP concentrations, while reservoirs had lower TP concentrations (Figure 5).Additionally, we observed that lakes and reservoirs with lower TP concentrations corresponded to lower Chl-a concentrations, as indicated by the median Chl-a concentration at the terminal nodes of the regression tree (5 mg L 1 for lakes and 2 mg L 1 for reservoirs).Chl-a concentration in lakes and reservoirs with low TP concentration was primarily influenced by TP concentration, as evidenced by the lack of further branching in the regression trees (Figure 5).These results are consistent with previous studies, including one that employed a process-based model, which demonstrated a positive linear relationship between TP and Chl-a concentration in a P-limited state, with the rate of linear increase determined by the C:P ratio of algal biomass (Olson & Jones, 2022).Similarly, a statistical analysis of 3,874 lakes in 47 countries worldwide revealed a sigmoidal relationship between TP and Chl-a, explaining 44% of the Chl-a variation, with a positive linear TP-Chl-a relationship at lower TP concentrations (0.004-0.23 mg L 1 ) (Quinlan et al., 2020).These findings collectively highlight the significant role of TP concentration in driving Chl-a variations in both lakes and reservoirs and provide valuable insights into the underlying mechanisms of Chl-a dynamics in aquatic ecosystems.
The regression tree analysis revealed distinct patterns in lakes and reservoirs with higher and lower TP concentrations (Figure 5).In lakes with higher TP concentrations, latitude, and longitude played significant roles in Chl-a variations, dividing the lakes into four spatial regions (Figure 5a, Figure S3 in Supporting Information S1).These regions displayed different mechanisms of influence on Chl-a concentration, with Region I (Yunnan-Guizhou Plateau) showing associations with TN and upstream watershed slope, Region II (Pearl River Basin) exhibiting complex effects, Region III (lower reaches of the Yangtze River and Yellow River basins, Haihe River Basin, and Northeast Plain) featuring diverse influences, and Region IV (middle and upper reaches of the Yangtze River and Yellow River basins and the Tibetan Plateau) demonstrating relatively low Chl-a concentrations (Figure 5a, and Figure S3 in Supporting Information S1).These findings corroborated previous research on freshwater lakes and reservoirs in China, which highlighted the significant influence of geographic location on the spatial variability of nutrients available to algae (Chl-a/TP and Chl-a/TN) (Huo et al., 2019).For reservoirs with higher TP concentrations, the slope of the upstream catchment emerged as the second most important predictor, dividing the reservoirs into high (≥6.41°)and lower (<6.41°)slope categories (Figure 5b).Furthermore, the average temperature in the months of April to September (T_AMJJAS) was the main influencing factor for reservoirs located in plains and low relief areas, while TP concentration remained the primary factor for reservoirs situated in mountainous areas (Figure 5b).The results suggest that the influence of slope on reservoirs could impact nutrient inputs due to differences in runoff and erosion rates, whereas plains are more affected by human activities, resulting in higher nutrient status and a greater influence of temperature on Chl-a (Huang et al., 2021;Paltsev & Creed, 2022;Zhou et al., 2022).
To further explore the different controlling factors on Chl-a in lakes and reservoirs with higher and lower TP concentrations, SEMs were developed (Figure 6).Four SEM models were constructed based on the MRT results, considering lakes with TP < 0.041 mg L 1 and ≥0.041 mg L 1 , and reservoirs with TP < 0.027 mg L 1 and ≥0.027 mg L 1 .The SEM analysis revealed the direct and indirect effects of nutrient concentrations, geographical characteristics, and meteorological factors on Chl-a concentration (Figure 6).In lakes and reservoirs with lower TP concentrations, TP concentration played a primary role in controlling Chl-a variations (Figures 6a  and 6c), consistent with the MRT results.Conversely, in lakes and reservoirs with higher TP concentrations, geographical and climatic factors exerted stronger influences on Chl-a than TP concentration (Figures 6b and 6d).Specifically, for lakes, geographical characteristics, extreme climate events, meteorological factors, and TN concentration had greater standardized coefficients for Chl-a than TP concentration (Figure 6b).Similar results were observed for reservoirs (Figure 6d).These findings highlight the importance of considering synergistic influences from geographical and climatic factors, along with nutrient concentrations, in understanding phytoplankton biomass dynamics when lakes and reservoirs are out of a P-limited state (Liang et al., 2020;Paltsev & Creed, 2022;Peng et al., 2021;Stockwell et al., 2020;Zhou et al., 2022).Extreme weather events, such as droughts, heat waves, and storms, can have transient and persistent effects on physical and biogeochemical processes in lakes and reservoirs, leading to changes in phytoplankton communities and dynamics (Frölicher & Laufkötter, 2018;Gobler, 2020;Michalak, 2016;Stockwell et al., 2020).

Implications for Lake and Reservoir Water Quality Management
Our comprehensive analysis of 255 lakes and 332 reservoirs in China revealed that the relationship between Chl-a and nutrients was not universally synchronous from 1980 to 2018.In particular, we observed that geographical characteristics, climatic factors, and TN concentration played a crucial role in influencing Chl-a concentrations when TP concentration exceeded a specific threshold (Figures 5 and 6).These findings challenge the traditional eutrophication management policies based on Liebig's law of minimum factors, which emphasize P control alone.Instead, our results support the idea that simultaneous control of both N and P is more effective for managing lake eutrophication (Conley et al., 2009;Ma et al., 2020;Paerl et al., 2016;Yu et al., 2019).
Geographical features exerted significant effects on Chl-a in lakes and reservoirs with higher TP concentrations (Figures 5 and 6).Given China's vast territory and complex geographical environment, there are substantial regional differences in the availability of algal nutrients in its lakes and reservoirs (Huo et al., 2019).Therefore, adopting a single nationwide set of nutrient standards and eutrophication control strategies may not be suitable for all regions.Instead, nutrient thresholds should be set based on nutrient ecoregions to achieve targeted water quality objectives and avoid exacerbating eutrophication in high-impact nutrient regions (Huo et al., 2019;Ma et al., 2023;Yu et al., 2019).
The influence of meteorological factors on the primary productivity of lakes and reservoirs was also notable (Figures 5 and 6).Climate change-induced alterations in weather patterns and variability have increased the frequency of extreme weather conditions, impacting water quantity and quality (Gobler, 2020;Merz et al., 2023;Michalak, 2016;O'Beirne et al., 2017;Stockwell et al., 2020).Consequently, it is imperative for China to adopt adaptive measures to protect water quality and lake and reservoir ecosystems.Implementing flexible pollution control strategies and advanced water treatment technologies to control nutrient loading are fundamental steps in nutrient management (Merz et al., 2023).Additionally, adopting a systematic governance approach that considers coordinated protection of water resources, water environment, water ecology, and water security is crucial for the holistic protection of lake ecosystems (Ma et al., 2023;Wang et al., 2022).This integrated approach will be instrumental in safeguarding China's water bodies in the face of potential climate change impacts and uncertainties.
Water Resources Research 10.1029/2023WR036587 ZHANG ET AL.

Conclusions
Our study contributes crucial insights into the changes in primary production and nutrient dynamics in lakes and reservoirs across China from 1980 to 2018 and the underlying influencing mechanisms.We achieved high accuracy in predicting Chl-a concentrations using a stacking model with an ensemble data set, with an NSE of 0.9658.Our findings reveal a concerning trend of increasing lakes and reservoirs at risk of algal blooms since 1980.At the individual lake and reservoir level, as well as at the regional level, we observed that the concentrations of TN and TP and their stoichiometric ratios did not always align with Chl-a concentrations.Notably, in lakes and reservoirs with lower TP concentrations (<0.041 mg L 1 for lakes and <0.027 mg L 1 for reservoirs), TP concentration emerged as the primary driver of Chl-a variations.However, when TP concentrations exceeded these thresholds, other factors such as geographical characteristics, climatic influences, and TN concentration played a critical role in shaping Chl-a concentrations.To effectively address and prevent eutrophication in lakes and reservoirs, it is essential to adopt dual control strategies for both nitrogen and phosphorus, tailored to specific nutrient ecoregions that account for regional variations in geographical features.Moreover, considering the uncertainties and risks posed by climate change, implementing flexible mitigation and adaptation strategies becomes paramount to safeguarding water quality and ecosystem health.Overall, our results provide essential information for informed and urgently needed decisions on managing water environments and addressing the challenges posed by eutrophication in China's lakes and reservoirs.

Figure 1 .
Figure 1.The methodological framework used in this study.
and phytoplankton biomass, changes in Chl-a, TN, and TP concentrations, and TN: TP in lakes and reservoirs in China over the past nearly 40 years were jointly analyzed.

Figure 2 .
Figure 2. Spatial variations in Chl-a concentrations from 1980 to 2018 in lakes and reservoirs of China based on predicated data using the stacking machine learning.

Figure 3 .
Figure3.Lake-specific and reservoir-specific rates of change for total nitrogen (TN), total phosphorous (TP), TN:TP and Chl-a were estimated using the mixed hierarchical linear model.Rates of change are only shown for lakes or reservoirs with trend estimates possessing 90% confidence intervals that did not overlap with zero.White circles represent the lakes or reservoirs that did not have significant trends.To compare the rate of change of different indicators, the scale values were unified between 5 and 5. Lakes and reservoirs with rates of change greater than 5 and 5 are represented by 5 and 5.

Figure 4 .
Figure 4.The relationship between estimated total nitrogen (TN), total phosphorus (TP) and Chl-a generated by the hierarchical linear model using the annual average value of nutrients and Chl-a concentration in lakes (n = 255) and reservoirs (n = 332) between 1980 and 2018.The vertical and horizontal lines in each panel represent zero change in the response variable on the x and y axes, respectively.The 45°line represents the 1:1 change ratio.In the contour plot (g, h), the change in Chl-a is shown by a color gradient.Numbers in the legend refer to the proportion (a-f) or number (g, h) of lakes or reservoirs in each symbol category or quadrant.

Figure 5 .
Figure 5. Multiple regression trees depicting the effect of nutrients, landscape and climate factors on phytoplankton biomass in 255 study lakes (a) and 332 reservoirs (b), using the monthly mean value of Chl-a concentration (μg L 1 ) as a proxy.N refers to the number of observations.

Figure 6 .
Figure 6.Structural equation model of phytoplankton biomass with nutrients, geography, soil properties, meteorological factors, and extreme climate in lakes (a, b) and reservoirs (c, d).The lake and reservoir data sets were classified according to TP concentration (0.041 mg L 1 and 0.027 mg L 1 , respectively) on the basis of the results from the regression tree analysis.Rectangular boxes indicate that the variable was be directly observable.The compound variable (shown as a circle) is a typical latent variable with error-free variance.Green arrows represent positive causal links, whereas red arrows indicate negative causal links.Dotted arrows indicate non-significant relationships, while solid arrows represent significant paths (P < 0.05).Standardized parameter estimates are shown next to each arrow.

Table 1
Estimated Parameters From the 3-Level Mixed Models, Where Intercepts Represent Concentration in 1980 and Percent Change per YearParameter ZHANG ET AL.