The Long-Term and Retention Impacts of the Intervention Policy for Cage Aquaculture on the Reservoir Water Qualities in Northern China

To ensure the safety of the water supply of the Panjiakou reservoir, in 2016, the Chinese central government comprehensively banned the fishing cage culture that had lasted for almost 30 years. However, the long-term effects and retention impacts of the government’s mandatory intervention on the reservoir water quality are unknown. To determine the reservoir water quality, we employed statistical methods along with the mathematical model to investigate the internal relationship since the construction of the reservoir. We applied seasonal trend decomposition using loess (STL) to explore the long-term and seasonality trend of monthly total nitrogen (TN) and total phosphorous (TP). To separate the impact of upstream water quality changes from cage culture on reservoir water quality, we employed generalized additive models (GAMs). We created a model, the LAKE2K model, to investigate the internal sources of the sediment that accumulated during the aquaculture period and its retardant effect. The results revealed that the concentration of upstream TN was more affected by non-point sources than by TP. The long-term policy of encouraging aquaculture has greatly contributed to the increase in the reservoir TP concentration rather than an increase in TN; the prohibition of cage aquaculture has resulted in a sharp drop in TP. After the ban, the sediment became the main source of TP. We suspect that the TP concentration of the reservoir and sediment will decrease gradually until a new equilibrium is reached within 10 years. This study offers lake managers an opportunity to increase their insight into the interaction of management measures with water quality and provides valuable information for the natural recovery of the eutrophic system.


Introduction
The tension of limited water supplies and increasing water demands is a common occurrence in many countries. To alleviate water shortage problems, inter-basin water transfer projects could be one of the effective engineering countermeasures by artificially re-allocating water resources [1,2].
Tianjin is an important industrial city in China, located next to Beijing, and is one of the most economically developed and densely populated areas; however, the city had suffered a severe water crisis in the 1980s-industrial discontinuation, domestic water limitation, water quality deterioration, and saltwater intrusion seriously affected the people's wellbeing, industrial production, and urban safety. To address the water shortage, China's central government launched the first cross-regional large-scale water diversion project of China in 1982 [3]. The Panjiakou (PJK) and Daheiting (DHT) reservoirs were built at the middle reach of the Luanhe River in Hebei Province and jointly operated for water storage and regulation ( Figure 1). Water from the cascade reservoirs subsequently flows into the Yuqiao reservoir in Tianjin [4,5]. The construction of this project significantly improved the quality of life for people in Tianjin. To address the poor quality-of-life of residents in the reservoir area, the local government began to encourage and support fish cage culture in the 1990s. After 2000, the area engaging in cage culture in the Panjiakou reservoir increased rapidly, exceeding 3300 hectares in 2014, with cages accounting for more than 6% of the water area ( Figure 2). As cage culture continuously expanded, a large amount of nutrient feed was put into the reservoir [6]. With the rapid development of cities upstream during the same period, the ecosystem of the reservoir experienced severe eutrophication, and the water quality has deteriorated significantly in the last decade, even affecting the normal water transfer plan and threatening the security of the water supply. Liu and Wang had calculated that nearly half of the total nitrogen (TN) annual loadings were from aquacultural pollution, where the contribution to the total phosphorous (TP) can reach up to 70% [7,8]. To protect water quality, the central government banned cage culture in 2016. As of 28 May 2017, more than 70,000 cages were removed. After the cages were cleared out, the water quality of the reservoir was comprehensively and rapidly improved [9]; however, whether the improvement of the water quality was caused by the removal of the cages is unclear, as upstream runoff may also have been responsible. This kind of government intervention, i.e., policies of long-term encouragement and then short-term mandatory prohibition, rarely occurs in the rest of the world; thus, it is necessary to explore the impact of these policies for cage aquaculture on the water quality of the reservoir. The availability of long-term monitoring data allows us to track the aquatic environment changes, and determine whether they are driven by natural and anthropogenic influences [10]. Deciphering these confounding factors often necessitates trend analysis [11,12]. Seasonal trend decomposition using loess (STL), a graphing technique, may be a powerful diagnostic tool [13,14]. It is flexible enough to reveal a much wider range of trend patterns and could provide visual results with seasonal interaction [13]. Furthermore, generalized additive models (GAMs) can produce regression functions that can be statistically described and compared [15,16]. It is constructed out of many smaller functions that can fit data with smooths or splines [17,18]. Both of these methods provide a flexible modeling framework to capture the complex, non-linear relationships among response and explanatory variables, which means that the potentially relevant covariates can be considered [12,19].
In addition, the Panjiakou reservoir has consistently accumulated fish feed since the 1980s; the nutrients derived from lake sediments may continue to mobilize long after external loadings have been reduced and retard ecosystem recovery [20,21]. Because of this, we must investigate the degree to which the retention of deposited fish feed affects water quality and how long it lasts [22,23].
There is a variety of techniques available for assessing internal nutrient loadings, such as in site measurements, laboratory experiments, mass balance methods, and others [24]. Field experiments provide snapshots of the situation, which is a highly dynamic process and which may create great uncertainty in reflecting the entire lake [23,25]. The results of the mass-balance method are based on an annual or seasonal scale [26]. In this regard, mathematical modeling, combining the diagenetic model with the ecosystem model, is a flexible and powerful tool to reproduce the dynamics of biogeochemical processes between sediment and bottom water [25,27]. LAKE2K is designed to simulate water quality trend considering geographic, physical, and biogeochemical conditions by Chapra and Martin [28,29]. This method contains a module to calculate oxygen and nutrient sediment-water fluxes, which can provide insights into dynamic sediment processes.
In this study, our aim was to investigate to what extent the long-term aquaculture and the sediment retention, after the dramatic management changes, could influence water quality. Our goals were to: (1) Identify the seasonal and long-term nutrient trends in the PJK reservoir and upstream inflow using the STL method since the construction of the reservoir. (2) Analyze the impact of longterm cage culture on the changes in reservoir water quality using the GAMs method. (3) Investigate how the accumulated fish feed in the sediment influences the water quality after the removal of the cage and to predict its retardant effect.

Study Area
The PJK and DHT reservoir system is located in northern Hebei province, China ( Figure 1). The PJK reservoir receives water from the upstream Luan River and two other small tributaries and delivers water to the DHT reservoir downstream. Then, the water is transferred to the Yuqiao reservoir in Tianjin and Tangshan city. These cascade reservoirs play a pivotal role in water supply, flood control, and power generation [30]. The length of the PJK dam is 1039 m with a maximum height of 107.5 m [4]. Its total storage capacity and utilizable capacity is 2.93 and 1.95 billion m 3 , respectively. The corresponding reservoir surface areas are 70 km 2 and 40 km 2 , respectively. In addition, the long term annual streamflow is about 2.45 billion m 3 and has transferred 19.22 billion m 3 water to Tianjin over 26 years [5,31]. The daily average temperature is between −14.7 and 31 °C over 30 years. The aquaculture area of the PJK reservoir has gradually increased from 1.54 hectares in 1988 and has entered a stage of rapid expansion after 2005. In 2014, the aquaculture area exceeded 3300 hectares, and cages accounted for more than 6% of the water area ( Figure 2).
The drainage area of Luan River above the PJK reservoir is 33,700 km 2 , which experiences a typical temperate continental monsoon climate. The annual precipitation and streamflow are unevenly distributed both spatially and temporally: the summertime precipitation accounts for more than 70%. The basin covers various landforms including forest (36.5%), grassland (34.1%), cropland (22.8%), etc. Although the proportion of urban land use is only 0.98%, there is a large city upstream and many factories along the Luan River, and because of this, the extensive use of agricultural fertilizers, the discharge of domestic sewage, and the production of industrial wastewater have a significant impact on water quality [3].

Data Sources
There are two monitoring stations in the study area: Wulongji (WLJ) station of Luan River at the mouth of the PJK reservoir and Panjiakou (PJK) station within the PJK reservoir. Periodic nutrient sampling began in the late 1980s, and regular monthly monitoring began in 2004 at the PJK station and 2015 at the WLJ station, respectively. The monthly data of TN and TP from 1988 to 2019 were used to reveal the long-term trends and explore possible causes. The hydrology data (inflow, outflow, water elevation, and water temperature), daily meteorology data (temperature, wind speed, and precipitation), and monthly water quality (Specific Conductance, Dissolve oxygen nitrogen, TN, and TP) in 2018 were used to establish the LAKE2K model. All raw data were provided by the Luan River Diversion Project Management Bureau. The data of the aquaculture area were obtained from the ChengDe Statistical Yearbook. Tables S1 and S2 show a detailed statistical breakdown of the dataset.

Statistical Methods
The STL method was applied to explore the long-term and seasonality trend of monthly TN and TP concentration at two stations. The STL method is a graphing technique based on various local smoothing techniques and has been widely applied [13]. Using this method, we can decompose a time series into three components as where Y, T, S, and R are the observed monthly value, trend component, the seasonal component, and the residual component, respectively; v is the measured data points. [14]. The decomposition is achieved through two loops. In the outer loop, robustness weights are assigned for the reduction or elimination of the effects of outliers. The inner loop iteratively updates the trend and seasonal terms using Loess (locally estimated scatterplot smoothing) smoothing [32]. The analysis was performed in R using the STL function. There are two smoothing parameters in the model, which represent the loess window degree of seasons and trend extraction, respectively, with window widths of 21 and 41, to best visually elucidate the trends [33]. We used the median polish method first to determine the missing monthly concentration values because the implementation of the STL method does not allow for missing data [34]. In the median polishing method, the trend component is estimated by the median value of the missing year and the seasonal component is estimated by the median of the missing month. The sum of these two items can estimate the missing value in the year and month. In our study, 29% and 10% of the monthly observations at WLJ and PJK stations required imputation, respectively. Figure S1 depicts the visualization of the long-term missing data using the "missingno" module in Python.
To separate the impact of upstream water quality changes from cage culture management measures on large interannual variation in reservoir water quality, the GAM method, proposed by Hastie and Tibshirani in 1990, was selected [35]. It is data-driven method rather than model-driven; that is, pre-specifying the form of the function is not necessary [17,36]. The model assumes that the mean of the dependent variable depends on an additive predictor through a non-linear link function. Suppose that Y is a response variable and X1, …, XP is a set of independent variables. The mathematical formulation can be expressed as: where E(Y) is the expected value of Y. The function g() allows a link between f(X1, …, Xp) and the expected value of Y, which amounts to an alternative for the response distribution besides just the normal distribution. μ is the model intercept. si(X), i = 1, …, p, are smooth functions wrapping the independent variable. In this study, we selected the gam() function from the "mgcv" package in R to fit a GAM. All water quality data were log-transferred before analysis to reduce the skewness of the variables. Six models were specified in the syntax as follows: gam2: gam(Panjiakou_P) ~ s(Area, k = 10), method = "REML")), gam3: gam(Panjiakou_P) ~ s(Wulongji_P, k = 10) + s(Area, k = 10), method = "REML")), gam4: gam(Panjiakou_N) ~ s(Wulongji_N, k = 10), method = "REML")), gam5: gam(Panjiakou_N) ~ s(Area, k = 10), method = "REML")), gam6: gam(Panjiakou_N) ~ s(Wulongji_N, k = 10) + s(Area, k = 10), method = "REML")), where "REML" is one of the smoothing parameter estimation methods instead of setting the smoothing parameter directly. k is the number of basis functions that make up a smooth function. Setting this value too low will prevent the model from being sufficiently wiggly. On the contrary, if it is high, the automatic smoothing parameter selection will prevent it from being too wiggly [37]. A good fit balances both over and under-fitting. "Area" is the area of cage culture within the PJK reservoir. Gam1 and gam4 used data from the WLJ station as predictors to fit the TN and TP concentration of the PJK station, separately. Gam2 and gam5 fit the TN and TP concentration of PJK station using the aquaculture area data as predictors, separately. Gam3 and gam6 considered these two straight-forward predictors concerning the relevant water-quality parameters. Before the exploration of model results, the gam.check() function was used to make sure that the models were well-fit. The comparison of model results can be evaluated by the goodness of fit statistics, including adjusted coefficient of determination (R 2 ), Akaike information criterion (AIC), and deviance explained ratio.

LAKE2K Model
In this study, LAKE2K (version1.4) was applied to model dynamics for several significant water quality variables of the PJK reservoir. The model was run within the Microsoft Windows environment: programmed in Windows macro language (Visual Basic for Applications, VBA) using Microsoft Excel as the graphical user interface [38]. The model runs on a daily scale and depicts the lake as a series of one-dimensional conical stratified lakes, consisting of three vertical layers (epilimnion, metalimnion, and hypolimnion)-see Figure S2. The change in the lake volume is computed through a dynamic water balance function: The heat balance is written for each vertical layer and ice thickness is computed at the air-water interface. Mass balance describes the sources and sinks of various model state variables (carbon, nitrogen, oxygen, phosphorus, silica, and phytoplankton and zooplankton biomass)-see Figure S3.
where ci is the concentration of layer i [mg/L or μg/L], cin is the concentration of the inflow [mg/L or μg/L], and S1 is sources and sinks of the constituent due to reactions and mass transfer mechanisms Furthermore, based on a model originally developed by Di Toro, the model can compute oxygen and nutrient sediment-water fluxes ( Figure S4). The approach divides the sediment into two layers: a thin (1 mm) surface aerobic layer underlain by a thicker (10 cm) lower anaerobic layer [29]. In the anaerobic layer, the diagenesis process happens: particulate organic matter settled from the overlying water is mineralized into the dissolved matter. These soluble reactive constituents are then transported to the aerobic layer where oxidization occurs. For a detailed description of the LAKE2K model, we refer readers to the Users' Manual and Documentation [29].
The model was simulated and calibrated for 2018. The input data include the daily hydrology data (Elevation-Area Curve, inflow, and outflow), daily meteorology (temperature, wind speed, and precipitation), and water quality (water temperature, DO, nitrogen, TN, and TP) of the WLJ station. Due to limited monitoring conditions, the input water quality data is on a monthly scale in some months. The low sampling frequency may result in uncertainty for the measured data. The interpolation calculation of the model on the input data has a significant impact on the results. Monthly water elevation, water temperature, DO, nitrogen, TN, and TP from the PJK station were used to calibrate the model parameters. Since water quality data were sampled on the water surface, the measured data were used to calibrate the model results in the epilimnion of the PJK reservoir. Figure S5 illustrates the inflow and outflow of the PJK reservoir in 2018. Since water quality components within the lake ecosystem are coupled and cross-pinning, model coefficients were adjusted artificially and iteratively until an adequate fit of measured versus simulated data was obtained. The root mean square error (RMSE) was used to quantitatively evaluate model performance for reproducing observed data [39]. It is calculated using the following equations:

RMSE
∑ P Q n where Qi and Pi are the observed and predicted values. n is the number of observed data. RMSE ranges from the optimal value of 0 to infinity. The closer the RMSE value is to 0, the better the performance of the model. Figure 3 illustrates the STL results of TN concentration at both stations. Due to cold weather and imperfect monitoring technology, the monitoring data in the winter (January, February, March, and December) of the previous few years were missing, especially at the WLJ station. Compared with monitored values, the imputed data have similar distributions; however, considering that there were more missing data at the Wulongji station, the credibility of the seasonal component trends in earlier years may be relatively low. Using the median polishing method to estimate missing data may artificially reduce the variance [40]. For the long-term trend component, both sites exhibited similar patterns (Figure 3c,d). After a period of modest change from 1988 to the late 1990s, TN concentration rapidly increased to the first peak around 2005, followed by a slight decline in the late 2000s, and then steadily increased to the historic peak. After that, the TN concentration dramatically decreased until 2017 and then a subsequent slight rebound occurred. The difference is that the TN concentration at the PJK station was lower than that at the WLJ station, and the upward and downward trends were more obvious. Furthermore, the fluctuation range of TN concentration at the WLJ station was significantly greater than that at the PJK station.

Long-Term and Seasonal Trend of Water Quality
For the seasonality patterns, TN concentration at two sites showed a distinct cycle. The seasonality of TN concentration at the PJK station was inconspicuous but gradually reinforced, indicating that the fluctuation of concentration within a year was gradually enhanced (left panel of Figure 3c). When the seasonal pattern was depicted by month, several features were highlighted (Figure 4e). TN concentration in summer and winter was slightly higher than for other months. Most months displayed fluctuations over the mean value, except June, July, and August, showing an increasing trend and December showing a decreasing trend throughout the records.
In contrast, the seasonality of TN at the WLJ station was steady and noticeable (left panel of Figure 3d), which followed a well-defined seasonal pattern: it rose steadily from January to June at its peak and decreased gradually until December (Figure 4f). This could be attributed to the influence of the regional monsoon climate in the basin and indicates that the upstream basin may be heavily polluted by non-point sources for TN. In addition, the long-term trend in each month was relatively stable over time. Figure 4 illustrates the STL results of TP concentration at both stations. Over the recorded period, the long-term smoothed TP concentration in both sites showed a different trend when compared to TN (Figure 4c,d). TP concentration rose continuously until it reached its peak in 2017, then sharply dropped to a relatively low value. Compared with the data of the WLJ station, the TP concentration at the PJK station had a larger magnitude of changes, which is the same as the difference in TN. These differences may be contributed to the combined effects of upstream inflow and aquaculture. The long-term seasonal tread of TP at PJK showed a similar trend to TN but had a higher range of variation (left panel of Figure 4c). In terms of monthly seasonal trends, TP showed a slightly pronounced seasonal cycle. There were high concentrations from April to August with an upward trend and low concentrations throughout October to May with a downward trend, which confirms that the fluctuation range has become larger (Figure 4e). This regular cycle differed from the TN and TP recorded at the WLJ station, which implies that the upstream basin has a small impact and might be associated with the feed from cage fishing and algae growth in spring and summer.
The seasonality of TP at the WLJ station was relatively stable but also unapparent (left panel of Figure 3d). Rather than the pronounced summer peak of TN, the monthly TP time series showed a distinct cycle: The long-term trend was relatively consistent across months with random fluctuation (Figure 4f). The TP in summer and winter was slightly higher than for other months, which suggests that the upstream basin may be heavily polluted by point sources for TP.
Overlaid with the original monthly time series, the STL estimated trend passes approximately through the middle of the observations and reproduced the main fluctuations in the data. Significant differences in the long-term trend and seasonality are observed between the two sites and the two water quality variables. Although the long-term trend of PJK tracked WLJ, the larger fluctuation at PJK may be contributed to the combined effects of upstream inflow and aquaculture. In addition, the turning point of TN was earlier than that of TP and has rebound in recent years. The long-term seasonality at PJK exhibited a reinforced tendency, whereas the seasonality of WLJ was relatively stable over time. The monthly seasonal pattern implies that the TN at the WLJ station was more influenced by the non-  [41]. It was concluded that nearly half of the TN annual loadings are from non-point sources, which mainly come from agriculture practices. TN loadings have a large seasonal variation, which coincides with our results. Generally, the best management practices (BMPs) are whatever practices are the most effective measure for controlling watershed pollution [42]. Qi et al. coupled the Revised Generalized Watershed Loading Function (RGWLF) model and the non-dominated sorting genetic algorithms II (NSGAII) optimization algorithm to identify the optimal spatial allocation of BMPs for dissolved nitrogen, which helps to realize cost-effective management [43].

GAM Results
According to the model check results, all the models converged after iterations. Table S4 summarizes the statistical results of model residuals. Figures S6-S11 plot the diagnostic four-panels for the residuals, including the Q-Q plot, residuals versus linear predictor, the histogram of residuals, and response against fitted values. All these outcomes demonstrate that the models performed satisfactorily overall, with a few exceptions for gam5. The gam5 diagnostic chart shows that the residuals are not completely randomly distributed, even with larger k values ( Figure S10), indicating that the aquaculture area cannot fully explain the changes in the TN of the PJK reservoir.
Fitting the data to each GAM model helps to detect the impact of each smooth term on the degree of fitting effect. Table 1 presents the results of the statistical comparison. Figure 5 illustrates the effect of the smooth functions of each variable. Not surprisingly, the increasing TP at the WLJ station tends to increase the TP at the PJK station. There is an obvious pattern in the aquaculture area series: a noticeable downward trend in the low values followed by a continued upward trend. Although the smooth function in both gam1 and gam2 was significant, gam2 performs well due to higher R 2 , deviance, and lower AIC. The aquaculture area could explain more than 25% of the variance than the upstream inflow. In addition, the model fitness seems to be improved by including both variables in gam3; however, the significance of s(log(Wulongji_P) is reduced at the 0.05 level and the value of effective degrees of freedom (EDF) is 1, which indicates a linear relationship. The results indicate that the TP at the PJK station is affected by the combined effects of upstream inflow and aquaculture, of which, aquaculture contributes more.  As for the models of TN at the WLJ station, the smooth function of s(log(Wulongji_N) has a positive and non-linear effect in gam4. The effect of aquaculture smooth function is similar to that of TP in gam5. According to the statistical results, smooth functions in both models are significant and the difference between them is noticeable. Gam4 performs well than gam5. The TN from upstream explained more than 96% variance and deviance of TN for the reservoir. Even including both smooth functions in gam6, the model performance only improved slightly. The s(log(area)) smooth term is no longer significant. These results imply that the TN from upstream determines the TN at the PJK station. According to the research results of Liu et al., the contribution ratios of cage fish culture to the TN and TP of the PJK reservoir are 53.1% and 70.9%, respectively. After the cages were cleared, non-point sources became the main source of pollution, and their contribution to TN and TP are 45% and 63%, respectively [8]. Aquaculture has a greater effect on TP than total nitrogen; thus, the banning of cages has significantly reduced the source of TP pollution and improved water quality [9].  Table 2 shows the RMSE results between the simulated and observed values. The RMSE values for all variables are less than 5%, indicating a very credible result when compared to previous studies [44,45]. Figure 6 illustrates the model calibration results of 2018 and shows that the model simulation outputs track the historical patterns and magnitude of water qualities. Although there are differences between the data and model predicted values, the results show that the model is ecologically interpretable. The water level rose in winter and spring for the water storage and dropped before the arrival of the flood season in summer. The water temperature data exhibit a parabolic shape, which ranged between 1 °C in winter and about 31 °C in summer. With the onset of spring, the DO steadily increased, which may be contributed to the available sunlight promoting the growth of phytoplankton. The variations of variables in the reservoir were consistent with the trends of upstream inflow. The water quality concentrations gradually declined before July, presumably influenced by the combined effects of settling, phytoplankton uptake, and lower upstream nutrients concentration; however, with the arrival of the summer rainy season, the abundant upstream pollutants and nutrients were brought into the reservoir, causing the TN and NH3 concentration to rise rapidly and then steadily decrease. Since phosphorus is less affected by non-point sources in the upstream basin, the increase in TP concentration in the reservoir was not as significant as the TN concentration.    [46]. The simulation results are in the same order of magnitude as the measured results. To predict the retardant effect of the sediment accumulated during the aquaculture period on the recovery of water quality after cage removal, the model was proceeded to run 10 years with the same boundary condition every year. Figure 8 shows the long-term trends of water quality and sediment and shows that the TP concentration trends steadily downward in the first 5 years and then stabilizes, whereas the TN gradually increases. This is consistent with our previous conclusions that the removal of fish in cages can significantly decrease the phosphorus concentration, and upstream water is still the main source of nitrogen. After external loads were substantially reduced, the exchange and content of sediment dropped rapidly in the first few years, because the P concentration in the reservoir decreased and sediment began to release P into the overlying water column. With the water mix and diversion, TP is generally drained out. Finally, a new equilibrium with a steady state between sediment and water TP concentration was achieved-TP concentration gradually approached 50 mg/L. The pattern of N behaves in the opposite manner due to the external N accumulation. Lewis et al. modeled the prediction of sediment response to changes in different fractions of downward P reductions and the results were 0.5 or 0.6 years for the fast fraction, where the release rate dropped rapidly, and 21 or 27 years for the slow fraction, where there was asymptotic progression toward steady-state release rates [27].

LAKE2K Model
Reducing the external nutrient load is the most crucial and direct way to restore the ecological quality of lakes and reservoirs besides bioremediation and physico-chemical methods. In some cases, obvious improvements were observed rapidly [47,48]. In other instances, internal loadings released from sediment delay the remediation [49][50][51]. A study examining 35 lakes revealed that in most lakes, a new equilibrium for TP was reached after 10-15 years [52]. Based on the actual conditions and model prediction results, the TP concentration in the PJK reservoir has been dramatically decreasing and the water quality and sediment of the PJK reservoir is expected to reach a new equilibrium within 10 years; however, due to the abrupt disappearance of filter-feeding fish, cyanobacteria have bloomed frequently in the past two years [8]. Such strong and sudden intervention management measures for cage fish culture may result in a drastic change in the structure and function of ecosystems, which could potentially cause irreversible damage to the state of the ecosystem [53]. In addition, the historical longterm use of antibiotics may have potential adverse effects on the ecosystem and public health [54].

Conclusions
In our study, we examined the long-term effects and retention impact of the government's mandatory intervention of cage aquaculture on the water quality of the PJK reservoir. The decadal and seasonal trend of monthly TN and TP concentration were graphically analyzed using the STL method. The GAM approach was applied to statistically identify the main factor affecting reservoir water quality from upstream water and cage aquaculture. The statistic results of model residuals and diagnostic four-panels demonstrate that the models performed satisfactorily. Additionally, we employed the LAKE2K model to reproduce the dynamics for several significant water quality variables of the PJK reservoir in 2018. Overall, the model was reliable for simulating water quantity and quality based on statistical analysis results. The following conclusions can be drawn: (1) Although the water quality changes in the reservoir were consistent with the upstream, there were still evident differences among the two sites and the two water quality variables, which implies the different influencing mechanism. It is concluded that the concentration of upstream TN was more affected by non-point source than TP.
(2) The inflow and aquaculture explained 97.2% variance of TN and 79.5% of TP for the reservoir, respectively. TN from upstream determines the TN at the PJK station. The results indicate that the TP at the PJK station is affected by the combined effects of upstream inflow and aquaculture, of which aquaculture contributes more.
(3) The long-term policy of encouraging the development of cage fish culture has greatly contributed to the increase in phosphorus concentration in the reservoir. The ban on cage aquaculture improved the water quality of the reservoir in a short time.
(4) The construction of the LAKE2K model indicated that the sediment is the source of the TP and its content decreased steadily after the removal of the fish cages. The results of the LAKE2K model also suggest that the TP concentration of the reservoir will decline gradually and a new equilibrium will be achieved within 10 years.
(5) Due to the greater influence of the upstream watershed, the concentration of TN shows an increasing trend, which requires more attention. Considering the limitations of the mathematical models, in-depth research and verification are needed in the future.

Supplementary Materials:
The following are available online at www.mdpi.com/2073-4441/12/12/3325/s1. Table  S1: Statistical information of the data used for trend analysis, Table S2: Statistical information of the data used for mathematical model, Table S3: Model state variables, Table S4: Statistical test for GAM in model residuals, Figure S1: Missing data visualizations (the more white lines in the matrix, the more missing values), Figure S2: LAKE2K water balance and vertical segmentation scheme, Figure S3: Model kinetics and mass transfer processes. Kinetic processes are hydrolysis (h), oxidation (x), nitrification (n), denitrification (dn), photosynthesis (p), respiration (r), death (d), grazing (g), and egestion (e). Mass transfer processes are reaeration (re), settling (s), sediment oxygen demand (sod), and sediment-water exchange (sw), Figure S4: Schematic of SOD-nutrient flux model of the sediments, Figure S5: The inflow and outflow of PJK reservoir in 2018, Figure S6: The residual diagnostic plot of Gam1, Figure S7: The residual diagnostic plot of Gam2, Figure S8: The residual diagnostic plot of Gam3, Figure S9: The residual diagnostic plot of Gam4, Figure S10: The residual diagnostic plot of Gam5, Figure S11: The residual diagnostic plot of Gam6.

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