Phosphorus Retention in Lakes: A Critical Reassessment of Hypotheses

Various hypotheses and models for phosphorus (P) retention in lakes are reviewed and 39 predictive models are assessed in three categories, namely mechanistic, semi-mechanistic, and strictly-empirical models. A large database consisting of 738 data points is gathered for the analyses. Assessing four pairs of competing hypotheses used in mechanistic models, we found that (i) simulating lakes as mixed-flow reactor is superior to plug-flow reactor hypothesis; (ii) modeling P loss as a second-order reaction outperforms the first-order reaction; (iii) P loss is better explained as a removal process throughout the lake volume than as a settling process across the sediments; and (iv) considering a fraction of P loading is associated with fast settling particles enhances lake total phosphorus (TP) predictions. Due to the systematic approach used for combining the hypotheses, some models are for the first time developed and assessed. For instance, the preeminent mechanistic model combines, for the first time, the second-order reaction hypothesis with the hypothesis that a specific proportion of P loading settles rapidly at the lake entrance. Results also showed that semi-mechanistic models outperform both mechanistic and strictly-empirical models since they take the form of a mechanistic model based on the physical representation of the lakes and utilize statistically acquired equations for unknown parameters. The best-fit model is a semi-mechanistic model that adopts the mixed-flow reactor hypothesis with a second-order volumetric reaction rate that is calculated as a non-linear function of inflow TP concentration, lake average depth, and water retention time. This model predicts 77.8% of the variability of log10-transformed lake TP concentration, which is 4.2% higher than the best mechanistic model and 0.8% higher than the best strictly-empirical model. The findings of this study not only shed light on the understanding of P retention in lakes but also can be useful for assessment of data-limited lakes and large-scale hydrological models to simulate the P cycle.


List of symbols
= Surface area of the lake ( 2 ) =Areal loading of TP (

96
A general TP mass balance model for the lakes, assuming that in the long-term the lake is estimated 97 as a Continuously Stirring Tank Reactor (CSTR), is as follows: Based on some previous models in the early 1960s and using the data of 8 Swiss lakes, 99 Vollenweider (1969) hypothesized that the loss of the TP from the lake water column to the 100 sediments is a linear function of the TP mass in water as follows: 105 1) can be rewritten as follows: By assuming that the mean water residence time ( ) in lakes is calculated as = ⁄ , 107 rearranging Eq. (3) takes the following form: where all the parameters except can be directly measured for a lake. Eq.
(3) assumes that there 109 are two outputs for the TP after entering the lake, i.e. it either is washed out of the lake or is retained 110 in the water column or is removed from lake volume via several reactions that are lumped and 111 simplified as a first-order reaction. However, other researchers (e.g., Chapra, 1975;Imboden, 112 1974; Lorenzen, 1973) treated the TP removal through the lake mainly as the sedimentation 113 process of P-containing particles with the settling velocity ( ) to the sediment surface (which is 114 assumed to be equal to lake surface area). In this approach Eqs.
directly measurable, Dillon and Rigler (1974) used the areal loading of TP ( , see Eq. 7) to 121 introduce the lake TP retention ( ) which is defined in Eq. (8).

122
= ⋅ (7) The input areal loading of TP is the sum of all the external loads of TP that enter the lake from 123 different sources and the output load is the output of TP loads through the lake outlet. Using this 124 approach, the loss term and the Vollenweider equation takes the following forms: literature as shown in next sections. 135 Prior research has interestingly enough suggested that empirical models of lake TP retention may 136 subsequently be explained with a mechanism. For instance, Jones and Bachman (1976) and is the TP concentration in cross-section . By simplifying and integrating Eq. (13), the 153 PFR lake model is as follows: where is the mean water retention time from lake entrance to cross-section . If is equal to 155 the length of the lake then = . By integration of Eq. (14), the mean is calculated as 156 follows: is linearly correlated to TP mass in the water column, which implies that the TP loss is the first-164 order reaction. This hypothesis was initially based on the data of four lakes in Vollenweider (1968).

165
Dillon (1974) theoretically investigated the use of a second-order reaction form. Walker (1985) 166 performed a more comprehensive study and investigated the case in which the loss term per unit 167 volume of the lake is a quadratic function of : is as follows: the mechanisms with lake TP retention (See Table 4 for their list), we decided to include them in 186 this study and assess the performance of all different types of models.

Materials and Methods
This section presents the materials, including the models and their classification criteria, and the 193 database of the lakes. The methods for fitting the models and their evaluation as well as the  Table   214 2), 13 semi-mechanistic (see Table 3), and 10 strictly-empirical models (see Table 4). Considering 215 that most of the semi-mechanistic and strictly-empirical models are non-linear, the refitting of the 216 models is conducted using the Genetic Algorithm heuristic search method in MATLAB 217 programming language (Appendix 1).
218 The effective loss rate 2 = The effective loss rate = The effective loss rate 2 =

230
The database used in this paper is a compilation of three data sets and has 738 observation data . Seventy-one lakes did not have data for and 149 lakes without 251 data for were also removed. Next, 5 lakes with a surface area greater than 10,000 km 2 (4 252 Laurentian Great Lakes and Lake Winnipeg in Canada) were excluded. Lake Tahoe in Nevada,

253
US, was also removed since its retention time ( = 700 ) is 11 times larger than the second 254 largest lake in the database ( = 60 for Lake Okanagan in British Columbia, Canada).

255
Considering that the net annual TP retention in lakes is assumed to be positive (i.e.

302
The overall best model(s) is(are) also highlighted in the last column.

306
The BIC estimate of mechanistic models (model #1-16 in Table 5) is used for the pairwise 307 comparison of the different hypotheses underlying the models. These comparisons include the 308 particle settling approach versus volumetric reaction approach; the hypothesis that lakes behave as  Based on the ∆ values (Fig. 5b), there is strong to very strong evidence that the mixed-flow 348 reactor hypothesis performs better than the plug-flow reactor hypothesis. In the literature, we found 349 only two studies that consider or compare these two hypotheses. Although Higgins and Kim (1981)   in Walker (1985). It is noteworthy to mention that the use of the second-order reaction models 376 does not add to the number of unknown parameters while increasing the prediction power. Another 377 difference between first-order and second-order models is that the second-order reaction model  i.e., models #11 and #15 can predict 16% and 23% respectively. However, models #11 and #15, 420 respectively perform about 2% and 8% better than models #9 and #13 in predicting the variation 421 of .
Using the ∆ < 2 criterion, the best intratype as well as the best overall models are chosen. a penalty is applied for the unknown parameters that result in negative simulated .

450
The overall comparison of the groups is also presented in Fig. 7. The semi-mechanistic models

483
The median and the Inter Quartile Range (IQR) of the relative errors are also presented in the corresponding panels.

485
The comparison of the mechanistic and semi-mechanistic models performance shows that the use is correlated to all three , and ̅ and as shown in Table 5, all four best performing models 498 are semi-mechanistic models whose TP removal rate is a function of these variables. The first-499 order and second-order volumetric loss rate of model #25 and #28 are as follows: order hypothesis is solved when using the dynamic loss term calculation in the semi-mechanistic 512 models, as the loss terms of first-order and second-order models in Fig. 9b are similar. However, 513 it is apparent that as the TP loss term increases (with the increase of lake TP concentration) the 514 behavior of the TP loss term in first-order and second-order models slightly differ. The first-order and a lower TP loss for lakes with higher .

524
The main objective of this paper was to assess four pairs of competing hypotheses that are 525 suggested for retention of TP in lakes using a large database. For this reason, 16 mechanistic 526 models are developed explicitly based on the physical representation of lakes. Specifically, this 527 research found that the best performing mechanistic model considers the lake as a mixed-flow 528 reactor where 30% of the input TP is rapidly settled in the entrance and the remaining participates 529 in a second-order reaction over the volume of the lake. It is worth highlighting that the -fraction 530 has been generally overlooked in previous studies and the combination of this hypothesis with 531 second-order reaction hypothesis and plug-flow models is for the first time conducted in this study.

532
Though the -fraction hypothesis is supported by the data, this fraction does not seem to be 533 constant for all lakes and this hypothesis overestimates TP retention for lakes with relatively short water retention time (e.g., < 1 month). Estimation of -fraction as a function of the lake and 535 river characteristics should be further investigated in the future. Using the lake and river 536 characteristics to calculate the unknown parameter of the mechanistic model results in the 537 development of a semi-mechanistic model, which is found to be the best performing type.

538
Modeling the TP removal as a second-order reaction outperformed the first-order reaction models 539 both in mechanistic and semi-mechanistic groups. The well-known strictly-empirical models not 540 only failed to perform better than the tested semi-mechanistic models but also they do not 541 necessarily provide any information about the retention mechanism. The results of this study 542 provide more insight into the P retention in lakes and can be used for large-scale hydrological 543 models to simulate P cycle and assessment of lakes eutrophication status. For finding the best models, the Bayesian Information Criterion (BIC) is used (Schwarz, 1978), 561 which take into account both the best fit and the number of calibrated parameters as follows: where the and are the indicator number of the model and, in this paper, is the model of lower 567 BIC estimate, i.e., the better model. By using the similarity to the likelihood ratio testing statistics, 568 Kass and Raftery (1995) have suggested the values in Table A.1 to be used for describing the 569 evidence against the model with higher BIC as a better model.