Watershed Response to Legacy Phosphorus and Best Management Practices in an Impacted Agricultural Watershed in Florida, U.S.A.

: Soil phosphorus (P) built up due to past management practices, legacy P, in the Lake Okeechobee Watershed (LOW) in south-central Florida, U.S.A., is often discussed as the root cause of lake eutrophication. Improvement of the lake’s water quality requires the identiﬁcation of critical P sources and quantifying their contributions. We performed a global sensitivity analysis of the Watershed Assessment Model (WAM), a common evaluation tool in LOW environmental planning, using the Morris method. A pre-calibrated WAM setup (Baseline) of the LOW sub-watershed, Taylor Creek Nubbin Slough (TCNS), was used as a test case. Eight scenarios were formulated to estimate the contributions of various P sources. The Morris analysis indicated that total phosphorus (TP) loads were highly sensitive to legacy P in improved pastures, the major land use covering 46.2% of TCNS. The scenario modeling revealed that legacy P, inorganic fertilizers, and other sources contribute 63%, 10%, and 32%, respectively, to the Baseline TP load of 111.3 metric tons/y to the lake. Improved pastures, dairies, citrus, and ﬁeld crops are the top TP load contributors. Our results have important implications for water quality improvement plans in the LOW and highlighted the need for accurate spatial mapping of legacy P and incorporation of such information in modeling efforts for watersheds demonstrating legacy P problems. for the LOW conditions needs assessment, followed by incorporation into agricultural and urban BMPs programs. The results of this study highlight the need for expediting projects and programs to help further mine


Introduction
The quantification of watershed nutrient loads in response to conservation practices, land use changes, and future climate projections is essential for developing sustainable watershed management plans [1,2]. Conservation programs that were implemented to reduce nutrient losses from non-point sources, widely referred to as the Best Management Practices (BMPs), have had only limited success in many North American, European, and Chinese watersheds [3][4][5][6][7]. Long-term accumulation of nutrients from the over-application of fertilizers or other prior land use practices (i.e., legacy nutrients) can lead to substantial lags in BMPs related water quality improvements and has been recognized as an important factor responsible for slower and lower than expected effectiveness of conservation programs in reducing nutrients [8][9][10][11][12][13].
The use of field and watershed-scale hydrologic and water quality simulation models for the development and assessment of environmental plans (e.g., Basin Management Action Plans-BMAP, Total Maximum Daily Loads-TMDL, Conservation Effects Assessment Project-CEAP) has become common practice [14][15][16][17][18]. While such models are extremely useful in estimating the potential impacts of BMPs, they do not capture legacy nutrient While P mass balance studies on LOW indicate that the net import of P was reduced by 25% since the late 1980s [35,50,51], studies of Lake Okeechobee inflow water quality found that TP loads going into the lake have not shown any decreasing trend since monitoring began in the 1970s, with an average annual value over 500 mt [52,53]. Improved pastures and dairies have consistently been reported as the top net P importers. A direct comparison of various mass balance components across LOW was not possible in the above studies due to differences in assumptions and methodologies used (e.g., sub-watersheds and basin delineation, study periods, mass balance components, etc.).
The Lake Okeechobee TP TMDL was adopted in 2001 and targeted the reduction in long-term LOW TP loads to 105 mt/y [54]. The deadline to achieve this TMDL was 2015, which was subsequently extended to 2035. Siders et al. [33], in a re-evaluation study of TMDL targets based on an enhanced water quality dataset for the lake, concluded that the original target of 40 µ g/L for the in-lake TP concentration is still appropriate. Lake Okeechobee's ability to assimilate phosphorus into its sediments has diminished over the last five decades, and the lake is estimated to turn from a P sink to a permanent source as early as 2035 [53]. This disconnect between environmental targets, actions, and current water quality warrants an investigation into the impacts of legacy P on TP loading. While the The Lake Okeechobee TP TMDL was adopted in 2001 and targeted the reduction in long-term LOW TP loads to 105 mt/y [54]. The deadline to achieve this TMDL was 2015, which was subsequently extended to 2035. Siders et al. [33], in a re-evaluation study of TMDL targets based on an enhanced water quality dataset for the lake, concluded that the original target of 40 µg/L for the in-lake TP concentration is still appropriate. Lake Okeechobee's ability to assimilate phosphorus into its sediments has diminished over the last five decades, and the lake is estimated to turn from a P sink to a permanent source as early as 2035 [53]. This disconnect between environmental targets, actions, and current water quality warrants an investigation into the impacts of legacy P on TP loading. While the Watershed Assessment Model (WAM) [55] was applied in formulating and updating the LOW BMAP, the focus has remained on assessing the effectiveness of source control programs. There has been a lack of attention towards formal WAM evaluation through uncertainty and sensitivity techniques, with only a single journal article focused on the routing module of WAM [56]. In order to improve the accuracy of the LOW BMAP, a better understanding of the legacy P representation and the relative importance of corresponding parameters is required. It is also important to understand the contributions of various nutrient sources for formulating and optimizing BMPs and other water quality improvement measures. The incorporation of legacy nutrient algorithms in watershed models is an active area of research, while only a few studies have apportioned nutrient loads to legacy and other sources with mechanistic watershed simulation models [19,57,58].
The aim of this study was to advance our knowledge of legacy P impacts using a hydrologic and water quality simulation model in a watershed with a recognized legacy nutrient problem, i.e., the Lake Okeechobee watershed. The specific objectives were to (a) perform a global sensitivity analysis (GSA) of WAM with the Morris method to assess the importance of legacy P and nutrient management parameters, and (b) quantify the relative contributions of various P sources (e.g., legacy, inorganic fertilizers, and others) and land uses to existing loads through hypothetical nutrient reduction scenarios. The Taylor Creek Nubbin Slough (TCNS) sub-watershed of LOW was used as a test case given its relative importance to the lake TP inflows, as noted in materials and methods.

Study Area
The TCNS sub-watershed is located on the northeastern shore of Lake Okeechobee and consists of five basins: S191 (488 km 2 ), S133 (104 km 2 ), S135 (72 km 2 ), S154 (129 km 2 ), and S154C (9 km 2 ). Together, these constitute only~7% of LOW ( Figure 1) but contribute disproportionately high TP loads to the lake. Currently, pastures make up the major portion of the sub-watershed area (~60%), followed by other agricultural activities (i.e., dairy, sugarcane, citrus, field crops, nurseries, and sod farms;~15.5%), natural land uses (wetlands and forests;~14%), and urban areas (~9.5%) ( Table S1). The four most abundant soil types include Immokalee, Myakka, Floridana, and Basinger, covering 35%, 21%, 13%, and 11% of the drainage area, respectively. All these soils are coarse-textured with high organic matter and are generally poorly drained [59]. A summary of the soil's physical and chemical properties is provided in the Supplementary Material (Table S2). Between the water years of 1991 to 2018, the TCNS had the highest TP flow weighted mean concentration (525 µg/L), the highest TP load (107 mt/y), and the highest TP load per unit area (1.33 kg/ha/y) among all the sub-watersheds discharging to Lake Okeechobee [60]. By considering these aspects, a variety of water quality improvement projects were implemented and planned for the TCNS, including dairy and pasture nutrient source control practices, wetland restoration, and stormwater treatment areas [56,61].

Watershed Assessment Model
The WAM is a physically-based, spatially distributed watershed-scale model that simulates the terrestrial part of hydrological and nutrient cycles [55]. WAM inputs include land use, soil, topography, stream networks, weather data (daily rainfall, monthly averageminimum and maximum temperature, solar radiation), and water control structures and their operations. The Basin Unique Cell Shell (BUCSHELL) sub-module divides the watershed area into square grid cells (typically 0.01 km 2 in size). Based on land use and soil type within a given model cell, BUCSHELL selects an appropriate field-scale hydrologic and water quality simulation model from a set of candidates: (a) the Everglades Agricultural Area Model (EAAMOD) [62][63][64] for agricultural and urban land uses on high water table soils, (b) the Special Case model [55] for all open water type land uses (e.g., wetlands, aquaculture, open pit mines), and (c) the Groundwater Loading Effect of Agricultural Management System (GLEAMS) [65] for all remaining land use and soil type combinations. The daily surface and sub-surface flows and nutrient loads generated by these field-scale models are adjusted to account for effects of impervious surfaces, wastewater systems, stormwater retention/detention features, and management practices that are not simulated mechanistically [66]. The Basin Land Area to Stream Routing (BLASROUTE) sub-module then routes the adjusted daily flows and nutrients from the source cells to the stream network and then to the watershed outlet. This routing is based on delayed unit hydrographs (defined separately for surface and sub-surface flows) for the source cell to the stream part and on a modified time-varying linear reservoir routing technique for the in-stream part [67]. Nutrient assimilation processes during source cell-to-stream and in-stream phases are based on generalized exponential decay equations that account for source-level concentrations, flow volumes, flow distances, and features such as wetlands and depressions encountered en route to the basin outlet [56].
where C = the outflow concentration (mg/L), Co = the inflow concentration (mg/L), Cb = the background concentration (mg/L), q = discharge (m 3 /s), a and b = the attenuation multiplier and exponent parameters for overland and stream water quality attenuation, d = the flow distance to a feature (m), ag = the attenuation multiplier for groundwater water quality, τ = the time step (typically 10 min), and R = the hydraulic radius (m). The parameters a, b, ag, and Cb are defined separately for each stream type, wetland type, and dryland (all non-wetland land uses) and for each water quality constituent. WAM is capable of representing diverse and unique hydrologic conditions in peninsular Florida, which are characterized by low gradients, abundant wetlands, high water tables, extensive drainage and irrigation canal networks, and complex water control operations. Most of the widely used watershed models, e.g., the Soil and Water Assessment Tool (SWAT) and the Hydrologic Simulation Program FORTRAN (HSPF), have limited capacities to routinely and precisely simulate one or more of these conditions without model improvements [68]. WAM was widely used in Florida for water quality assessments in research and restoration planning [69][70][71][72][73].

Taylor Creek Nubbin Slough Modeling
In this study, we used the TCNS WAM setups obtained from the Florida Department of Environmental Protection (FDEP). They were developed as a part of the Lake Okeechobee BMAP, the latest 5-year update, to more accurately depict existing conditions and potential impacts of restoration projects [61]. The first WAM setup depicts the existing conditions, Baseline from here on, which has undergone a multi-gage calibration and validation (2003-2013) [63] (Table S3). It incorporated currently implemented source control programs and their spatial extent. The performance of the Baseline WAM in simulating the TCNS wide flows and TP loads over the combined calibration-validation period is shown in Figure 2. The Nash-Sutcliffe Efficiency values for flow and TP loads were 0.85 and 0.74, respectively, while the percent biases were +4.17% (underprediction) and −2.92% (overprediction), respectively. This implies that the Baseline model's performance was satisfactory. Model performance for individual basins is presented in the Supplementary Material (Table S3 and Figure S1). The second WAM setup obtained from the FDEP assumed the implementation of advanced fertilizer management practices on all agricultural and urban lands. The management practices included fertilization rates and/or application timing recommendations made in the most recent Florida Department of Agricultural and Consumer Services (FDACS) BMP manuals [74]. This alternative setup is called FBMP in this article. To assess the relative impacts of various nutrient reduction projects, SWET [74] compared the 21-year (1993-2013) simulation results of restoration alternatives with the Baseline; we adopted the same approach in this study. The FBMP setup was used only indirectly to assign parameter ranges in the sensitivity analysis and to create some of the new scenarios for the segregation of P source contributions. Under the current spatial distributions of land uses and soils in the TCNS, WAM simulated 86.4% of the source cells in this sub-watershed with the EAAMOD, while the Special Case and GLEAMS were selected for only 10.6% and 3% of the TCNS domain, respectively. Clearly, the processes represented by the EAAMOD and their respective parameterizations are important for the accuracy of WAM's application to the formulation of nutrient source control plans for the TCNS. Under the high water table conditions in South Florida, the majority of flow and P transport occurs as surface and shallow subsurface flows. Therefore, P levels in the uppermost soil layers are the most influential in P transport [39]. A number of field-scale hydrologic and P transport models were developed in the 1990s to simulate these unique conditions, amongst which the EAAMOD was identified as the most suitable [75]. The EAAMOD underwent updates to enhance its usage from agricultural parcels on Histosol soils to all high water table soils with agricultural and urban land uses [63,76].
The EAAMOD includes a simple conceptual phosphorus module that accounts for various P pools and transport processes, including P in rainfall, irrigation and fertilization, removal by crop uptake, vertical movement with water (leaching and ET), mineralization of organic P depending upon the saturation condition, horizontal P movement with flow, and P partitioning between soluble, adsorbed, and suspended phases (when flooded). The EAAMOD's P module works on the same spatial discretization framework as flow. The estimation of the P mass balance is closely related to the water inflows and outflows and to the hydraulic properties and variables. The aerobic/anaerobic mineralization and the P adsorption/desorption processes depicted in the EAAMOD are important to P dynamics. Additional technical details on the EAAMOD flow and P modules are provided in the Supplementary Materials. A unique feature of the EAAMOD is that it explicitly accounts for the effects of TP accumulation from prior land use/management practices based on userspecified initial P application to the top soil layer [67]. The EAAMOD tracks changes in P forms and concentrations in the soil profile over the simulation period. The Baseline WAM setup determined legacy P concentrations for various land uses based on a single typical value for each land use obtained from the literature [40,63]. In a restoration assessment with WAM, the legacy P parameter values are reduced to account for the anticipated reduction in legacy P under the implementation of source control programs. This reduction is typically determined from the net P export for given land use [73,77].

Global Sensitivity Analysis
The Morris method [78], which is a GSA technique, was employed to assess the relative importance of phosphorus sources, legacy P, and inorganic fertilizer P on TP loads outflowing from the TCNS. It is one of the most extensively used sensitivity analysis techniques for computationally intensive and complex large-scale system simulation models, including watershed-scale hydrologic and water quality modeling tools, due to its global nature and ability to accurately and efficiently identify parameter importance with relatively very few model runs [79][80][81][82]. The Morris method is based on the calculation of multiple local derivatives (Equation (4)) for each of the k model parameters of interest. These derivatives are popularly called elementary effects (EE) in the literature on the Morris method. The parameter sample consists of r independent trajectories. Each trajectory consists of (k + 1) points traced within the parameter hyperspace such that one EE is evaluated per parameter from it. EE distributions are further analyzed to calculate three sensitivity measures for each parameter, µ* (Equation (5)), µ (Equation (6)), and σ (Equation (7)).
where y is the model output for the corresponding parameter set, x i is the ith input parameter, EE i is the elementary effect for the ith input parameter from a given trajectory, ∆ is the jump size by which x i is changed along the ith parameter axis, and r is the number of parameter trajectories. Parameters are plotted in µ* − σ and µ − σ planes to identify their relative importance and nature of impact (monotonic vs. non-monotonic; interactions vs. additive effects) on the output/s of interest [56,83]. µ* was shown to be an indicator of the overall impact of a parameter on model output, while σ indicates the interactions [84]. Parameters that are well separated from the origin of µ* − σ plane are important, while ones close to the origin are non-influential and can be fixed for further model analysis such as calibration.
The number of model runs required for the Morris method is estimated as r(k + 1). There is a range of sampling techniques associated with the Morris method that differ in principle [85]. A p-level fixed grid approach, such that ∆ = p/[2(p − 1)], is the most commonly used [86]. While r = 10 to 20 and p = 4 are recommended [87,88], Chitale et al. [89] showed that selecting an r that is a multiple of p provides more accurate results at a lower r.
By considering the purpose of this study (i.e., identification of the relative importance of legacy P and fertilizer P parameterization), the Morris method analysis was restricted to the EAAMOD simulated land uses that covered at least 0.5% of the TCNS (Table 1). Overall, 15 parameters (8 legacy P and 7 inorganic fertilizer P parameters) representing 11 land uses were selected. The enhanced Sampling for Uniformity algorithm [89,90] was used with p = 4 to generate 16 trajectories to form 16 × (15 + 1) = 256 parameter samples. All parameters were assumed to have uniform distributions with values representing reductions (expressed as fractions) from the Baseline TCNS WAM setup. While parameter probability distributions can have impacts on sensitivity analysis results, assuming uniform distributions is recommended when only the basic information (uncertainty range) is available [88,91]. The upper limits of the parameter uncertainty ranges were assigned based on the values for the corresponding parameters in the FBMP setup [74] and past WAM applications for restoration assessments [70,73].

Scenario Analysis for Source Type and Land Use Contributions
The relative contributions of various P sources such as legacy P, inorganic fertilizers, and all remaining P sources (mainly organic fertilizers/animal manure, etc.) and important land uses were quantified in the scenario analysis. In addition to the Baseline, eight alternative scenarios were formulated by separately reducing the legacy P and fertilizer P parameter values from the Baseline to those in the FBMP setup and finally to 0 ( Table 2). All agricultural and urban land uses where legacy P and/or fertilizer P can have an impact were incorporated into the scenario analysis, unlike the Morris method GSA that was limited to the EAAMOD simulated dominant land uses. Thus, the scenario analysis included both the EAAMOD and GLEAMS. A summary of legacy P and inorganic fertilizer P parameter values by land use are provided in the Supplementary Material (Table S4). Organic fertilizer/manure application rates and rainwater P concentrations were not changed in any of these scenarios. The above scenarios were focused on changes in source-level P losses. Physically, changes in source-level nutrient loads would impact nutrient assimilation in overland and in-stream flow phases. In WAM, nutrient assimilation processes are represented by Equations (1) and (3) and any changes in assimilation capacity would be rendered by changes in the background concentration (Cb), attenuation multiplier (a), and attenuation exponent (b) parameters. However, in WAM, all these parameters are user-defined. Khare et al. [56] showed that in WAM's routing sub-module, P dynamics and loads at the watershed outlet are most sensitive to attenuation multipliers for streams and dominant land uses. On the other hand, Cb parameters are estimated based on observations and would be reduced as a result of nutrient flushing. These parameters warrant consideration while quantifying the potential changes in TP loads at the watershed outlet under source-level reductions. A precise estimation of the impact of source-level phosphorus reductions on these parameters is uncertain. Therefore, eight alternative scenarios were further quantified with three levels of changes in Cb and a parameters-low, medium, and high, resulting in nine variations for each. These levels for Cb (a reduction from the Baseline by 0%, 10%, and 20%) and for a (an increase from the Baseline by 0%, 50%, and 100%) were chosen to demonstrate assimilation scenarios and were based on the authors' best professional judgment. Baseline values for these parameters are provided in the Supplementary Material (Table S5).
Median values of TP loads for the eight scenarios obtained from the above nine variations in assimilation capacity were used to calculate the contribution of legacy P, inorganic fertilizer P, and the remaining P sources (e.g., organic fertilizer/manure, rainfall, etc.). For example, the difference between the Baseline and ZeroLP, BaseFP scenario (i.e., zero legacy P and Baseline inorganic fertilizer P) would be an estimate of the contribution of legacy P from major agricultural land uses to the Baseline TP loads. Similarly, the difference between the Baseline and BaseLP, ZeroFP scenario (i.e., Baseline legacy P and zero inorganic fertilizer P) would represent the potential impact of inorganic fertilizer applications. The ZeroLP, ZeroFP scenario (zero legacy P and zero inorganic fertilizer P) would represent the combined contribution of organic fertilizers, manure applications, P in rainfall, and P losses from wetlands and other natural land uses. The RedLP, RedFP scenario (reduced legacy P and reduced inorganic fertilizer P) is equivalent to the FBMP abatement alternative in the LO BMAP evaluation and represents the plausible future TP loads from the TCNS when advanced nutrient management BMPs are effectively implemented. The BaseLP, RedFP (i.e., Baseline legacy P and reduced inorganic fertilizer P) and RedLP, BaseFP (i.e., reduced legacy P and Baseline inorganic fertilizer P) represent the contributions of reductions in inorganic fertilization rates and potential reductions in legacy P corresponding to these fertilizer reductions, respectively. The relative importance of major land uses was assessed by summarizing the source-level unit area TP loads (UATPL) and total source-level TP load contributions under the Baseline and eight scenarios.

Global Sensitivity Analysis
The sensitivity of average annual TP loading to Lake Okeechobee was illustrated by plotting legacy and fertilizer phosphorus related model parameters (Table 1) in µ* − σ and µ − σ planes (Figure 3a,b). The parameters' importance ranks were evaluated based on the µ* sensitivity statistic. The average annual TP load (sum of TP loads from all TCNS outlets) was most sensitive to the improved pasture legacy P parameter, indicating its critical role in P dynamics for the TCNS, followed by citrus fertilizer phosphorus. The field crop fertilizer phosphorus was ranked third but was not as separated from the origins of the µ* − σ and µ − σ plots compared to the improved pasture legacy P and citrus fertilizer P. The remaining parameters were relatively close to the origin on both µ* − σ and µ − σ plots, indicating they had little influence on TP loads. The dotted lines in these plots, µ* = σ in µ* − σ and µ = ±2 SEM = ±2 σ/sqrt(r) in µ − σ are thresholds to identify the degree of parameter interactions. Parameters that are located above the µ*= σ line on the µ* − σ plane or inside the "V" formed by µ = ±2 SEM on the µ − σ plane exhibit stronger interactions compared to their direct effects on the model output. Figure 3a,b indicate that σ values for all parameters were small, and they are located well below their respective threshold lines in both plots, indicating minimal to no parameter interactions. Another effect identified from the Morris GSA plots is whether the relationship between the model outputs and the parameters is monotonic or non-monotonic. This is verified by assessing µ and µ*. The relationship is monotonic if µ* = abs(µ). For the TCNS, the average annual TP load was found to behave monotonically with respect to all parameters. Moreover, for all parameters, the µ values were negative, indicating that the TCNS TP loads decreased with larger reductions in legacy P and fertilizer P. lets) was most sensitive to the improved pasture legacy P parameter, indicating its critical role in P dynamics for the TCNS, followed by citrus fertilizer phosphorus. The field crop fertilizer phosphorus was ranked third but was not as separated from the origins of the µ* − σ and µ − σ plots compared to the improved pasture legacy P and citrus fertilizer P. The remaining parameters were relatively close to the origin on both µ* − σ and µ − σ plots, indicating they had little influence on TP loads. The dotted lines in these plots, µ* = σ in µ* − σ and µ = ±2 SEM = ±2 σ/sqrt(r) in µ − σ are thresholds to identify the degree of parameter interactions. Parameters that are located above the µ* = σ line on the µ* − σ plane or inside the "V" formed by µ = ±2 SEM on the µ − σ plane exhibit stronger interactions compared to their direct effects on the model output. Figure 3a,b indicate that σ values for all parameters were small, and they are located well below their respective threshold lines in both plots, indicating minimal to no parameter interactions. Another effect identified from the Morris GSA plots is whether the relationship between the model outputs and the parameters is monotonic or non-monotonic. This is verified by assessing µ and µ* .The relationship is monotonic if µ* = abs(µ). For the TCNS, the average annual TP load was found to behave monotonically with respect to all parameters. Moreover, for all parameters, the µ values were negative, indicating that the TCNS TP loads decreased with larger reductions in legacy P and fertilizer P. In order to assess the impact of the spatial distribution of land use within the TCNS, the Morris method analysis was conducted for the average annual TP loads for the individual TCNS basins (i.e., at S191, S154, S133, S135, and S154C). Results are presented in Figure 4 a-j. While the model behaviors (monotonic response and low parameter interactions) observed for entire TCNS average annual TP loads were also evident at the basin level average annual TP loads, the parameter importance ranks differed. In the case of S191, the parameter importance ranks were similar to those of TCNS, with the top four ranked parameters being identical (Table 3). This is because S191 is the largest basin by area (~61%) and the largest contributor of TP loads (over 70%) among the TCNS basins ( Figure S2). In S154, the improved pasture legacy P was the most influential parameter, while the abandoned dairy legacy P and nurseries fertilizer P were ranked second and third, respectively. However, they were relatively closer to the origin of the Morris plots, indicating a relatively low influence on model output. Sod fertilizer P was the top-ranked in S133, followed by improved pasture legacy P and field crop fertilizer P. S135 exhibited four influential parameters that included sugarcane fertilizer P, nurseries fertilizer P, In order to assess the impact of the spatial distribution of land use within the TCNS, the Morris method analysis was conducted for the average annual TP loads for the individual TCNS basins (i.e., at S191, S154, S133, S135, and S154C). Results are presented in Figure 4a-j. While the model behaviors (monotonic response and low parameter interactions) observed for entire TCNS average annual TP loads were also evident at the basin level average annual TP loads, the parameter importance ranks differed. In the case of S191, the parameter importance ranks were similar to those of TCNS, with the top four ranked parameters being identical (Table 3). This is because S191 is the largest basin by area (~61%) and the largest contributor of TP loads (over 70%) among the TCNS basins ( Figure S2). In S154, the improved pasture legacy P was the most influential parameter, while the abandoned dairy legacy P and nurseries fertilizer P were ranked second and third, respectively. However, they were relatively closer to the origin of the Morris plots, indicating a relatively low influence on model output. Sod fertilizer P was the top-ranked in S133, followed by improved pasture legacy P and field crop fertilizer P. S135 exhibited four influential parameters that included sugarcane fertilizer P, nurseries fertilizer P, citrus fertilizer P, and improved pasture legacy P, in reducing order of importance. S154C had improved pasture as the only land use type, and the corresponding parameter was found to be the only influential parameter. Overall, this variability in parameter ranking was closely related to the relative abundance of land uses in the respective basins (Supplementary Material, Table S1).

Parameter
TCNS S191 S154 S133 S135 S154C  Figure 5 summarizes the TCNS TP loads under the Baseline and eight alternative scenarios. The Baseline scenario represents the calibrated-validated model. For simplicity, the contributions of different P sources were calculated based on their median values for each scenario. Results for legacy P, inorganic fertilizer P, and all other sources for the TCNS and individual basins are summarized in Figure 6. The contributions across the TCNS were estimated to be: legacy P-70.2 mt/y (67.7 mt/y to 70.10 mt/y), inorganic fertilizer P-17.53 mt/y (11.2 mt/y to 23.4 mt/y), and all other P sources-35.1 mt/y (33.2 mt/y to 37.2 mt/y), representing 63% (61% to 65%), 16% (10% to 21%), and 32% (30% to 33%) of the Baseline loads, respectively. It is worth noting that the average 17.53 mt/y (16%) contribution of inorganic fertilizer P is an artifact of unrealistically large changes in the Cb and a parameters for the corresponding source reductions. We expect inorganic fertilizer contributions to be close to 11.2 mt/y (10%), which corresponds to the Baseline values of Cb and a. Since inorganic fertilizer P-based TP load reductions are small, changes in their assimilation capacity would also be small. Using the same logic, the mere reduction in inorganic fertilizer P rates under advanced nutrient BMPs (i.e., FBMP) would reduce the TCNS TP loads by 8 mt/y while corresponding changes in legacy P under FBMP would be responsible for 4.4 mt/y. This implies that a considerable amount (approximately one-third) of the TP load reductions anticipated under FBMP (RedLP, RedFP) are attributable to legacy P uptake estimates.
The sum of the contributions from the three sources does not add up to the Baseline loads. First, the parameters related to nutrient losses from wetlands and other land uses corresponding to the Special Case field-scale model were not varied in this assessment. That means that source-level loading from 10.6% of the TCNS remained the same across all eight alternative scenarios. Second, the overland and in-stream nutrient assimilation conceptual model in WAM modulates the source-level nutrient concentration to the watershed outlet. As explained earlier, we created secondary scenarios by changing the Cb and a parameters. These are known to interact in a complex manner [56] and can contribute to source apportioned estimates that do not sum to Baseline loads. A similar observation was made by Kast et al. [58].
In the case of the S191 basin, the percentage source contributions were similar to the TCNS, as expected. Similar to the Morris method analysis results, the other TCNS basins exhibited relative source contribution estimates different from the overall TCNS due to different land use distributions. For example, in the S135 basin, sugarcane is the most abundant land use (~27.5%) along with relatively higher coverage of tree nurseries (~10%) and citrus (~3.5%), resulting in a higher impact of inorganic P fertilizers.
Cb and a. Since inorganic fertilizer P-based TP load reductions are small, changes in their assimilation capacity would also be small. Using the same logic, the mere reduction in inorganic fertilizer P rates under advanced nutrient BMPs (i.e., FBMP) would reduce the TCNS TP loads by 8 mt/y while corresponding changes in legacy P under FBMP would be responsible for 4.4 mt/y. This implies that a considerable amount (approximately onethird) of the TP load reductions anticipated under FBMP (RedLP, RedFP) are attributable to legacy P uptake estimates. Figure 5. Average annual Total Phosphorus loads from TCNS to Lake Okeechobee under the Baseline conditions and eight hypothetical scenarios of reduced legacy P and/or inorganic fertilizer applications ( Table 2). In each boxplot, the box represents the interquartile range, the whiskers show the 95% confidence limits, the middle horizontal line shows the median value, and the cross sign represents the arithmetic average.  (Table 2). In each boxplot, the box represents the interquartile range, the whiskers show the 95% confidence limits, the middle horizontal line shows the median value, and the cross sign represents the arithmetic average. The above findings of the relative impacts of P sources were qualitatively supported by the source-level UATPL maps (Figure 7), where areas in shades of red and green indicate relatively large and small source-level unit area loads. The reduction in legacy P has more impact on UATPL (moving from left to right) compared to the reduction in inorganic fertilization (moving from top to bottom). Table 4 summarizes the TCNS average annual source-level UATPL for major land use categories along with their contributions to the overall source-level TP loads. Dairy-related land uses (active dairies, outer pasture, and intensive pastures) had the top UATPLs with TCNS average values ranging from 4.71 kg/ha/y to 11 kg/ha/y across all nine scenarios. Their overall contribution to source-level total TP loads was the second highest (16.5% to 29.7%). Improved pastures contributed 30% to 54.5% of TP loads to these scenarios, the highest of all land uses, while the corresponding UATPL contributions were relatively low (0.43 kg/ha/y to 1.82 kg/ha/y). Other land uses that are important based on the UATPL include abandoned dairies, row crops, spray fields, and other agricultural activities (e.g., various types of animal feeding operations) with values in excess of 4 kg/ha/y. Among these, abandoned dairies contributed over 10% of source-level loads for six of the nine scenarios. Overall, legacy P was found to contribute 87.5 mt/y of the TP load at the source-level, of which improved pastures, abandoned dairies, and active dairy activities were responsible for 49.5 mt/y, 14.5 mt/y, and 11.6 mt/y, respectively. Inorganic fertilizers contributed 18.6 mt/y of source-level TP loads, of which a major portion (6.7 mt/y) was associated with intensive and outer pastures in the dairy land use class. Citrus (2.8 mt/y) and field crops (1.8 mt/y) were the second and third highest TP load contributors based on their inorganic fertilizer sources. Improved pastures and dairies constituted approximately 69% of the 41.8 mt/y source-level loads that originated from all other (mainly manure) P sources. The source-level contributions quantified from the scenario analysis are consistent with the results of the Morris method GSA, with the exception of active dairy land uses. Active dairies were not included in the Morris method GSA due to their small coverage and/or no anticipated changes in legacy P or fertilizer P values as depicted in the FBMP setup. The sum of the contributions from the three sources does not add up to the Base loads. First, the parameters related to nutrient losses from wetlands and other land corresponding to the Special Case field-scale model were not varied in this assessm That means that source-level loading from 10.6% of the TCNS remained the same ac all eight alternative scenarios. Second, the overland and in-stream nutrient assimila conceptual model in WAM modulates the source-level nutrient concentration to the tershed outlet. As explained earlier, we created secondary scenarios by changing th and a parameters. These are known to interact in a complex manner [56] and can con ute to source apportioned estimates that do not sum to Baseline loads. A similar obse tion was made by Kast et al. [58].

Hydrologic and Water Quality Simulation Models and Legacy P
Water quality impairment and other environmental consequences originating from the intensive use of fertilizers in agriculture for increased food and energy production have intensified globally in recent years [6,92]. High fertilizer applications have also resulted in substantial increases in soil nutrient stores. In many watersheds in the U.S.A., these legacy build-ups are the dominant source of current riverine nutrient exports [12,93]. The watershed-scale quantification of nutrient dynamics, including legacy nutrients, is an important facet for developing effective nutrient management strategies for improving water quality, but spatial-temporal heterogeneities make it difficult to accomplish this through direct monitoring efforts. Both statistical and mechanistic models are used for water quality assessments but need improvements to reduce uncertainties and accurately and reliably depict residence times and biogeochemical processes related to legacy nutrients [11,19,93]. For example, some of the drawbacks include (a) inaccurate assumptions about the export times in the case of statistical models, (b) a lumped parameterization to represent legacy nutrients effects, and (c) a lack of efficient model validation at the field-scale for mechanistic models. However, these models can provide valuable insights and are useful research and decision support tools.
A comprehensive summary of legacy nutrient modeling by Chen et al. [19] and our own literature search indicated that there are only a few modeling studies that have assessed legacy P effects with mechanistic watershed-scale hydrologic and water quality models. In the Maumee River Watershed, a large Lake Erie tributary, legacy P impacts were evaluated using SWAT with a scenario analysis approach [58,94]. While Muenich et al. [94] focused on calculating how long legacy P alone would be able to sustain crop growth, Kast et al. [58] estimated that legacy P contributed 44% to the outflow TP loads. Han et al. [57] found that SWAT outputs are sensitive to parameters relevant to legacy P dynamics in the agriculture-dominated Dengsha River Watershed in China. In the Yahara Lakes Watershed in Wisconsin, another agriculture-dominated basin, phosphorus losses were found to be sensitive to legacy P with an estimated load contribution of 35% [95]. An integrated modeling approach linking the Agro-IBIS (a processes-based ecosystem model), THMB (a flow, sediment, and nutrient routing model), and the Yahara Water Quality Model (a lake water quality model) was used in that study. A common link in the above studies was a lumped parameterization to portray legacy P effects and a model sensitivity assessment through scenario analysis. Our study performed a quantitative assessment to answer legacy P-related questions with WAM, a process-based watershed-scale model, in a region that has exhibited high agricultural legacy P accumulation in the U.S.A. [25]. We applied a more advanced global sensitivity analysis technique and a scenario analysis to explore the relative importance of legacy P and fertilization, specifically considering land uses in both analyses. In addition to WAM's suitability for simulating typical hydrologic conditions in Florida, its explicit representation of legacy P made this analysis possible.
Legacy P is essentially a result of the application of phosphorus in excess of plant requirements due to the use of fertilizer and manure that builds up in soils over decades [95,96]. Consequently, the legacy P at a given location is closely tied to the soil properties, land use/land cover, management practices, and temporal changes that occur in the latter two factors [97][98][99]. Typical approaches used in estimating terrestrial legacy P for informing quantitative models include mass balances, chemical analyses, and the phosphorus cycle [100]. Lou et al. [100] pointed out the drawbacks of the above methods in accurately identifying spatial variability in legacy P, demonstrating the need to combine traditional methods with remote sensing and advanced geostatistical techniques. In WAM (EAAMOD), the legacy P parameter values are based on land use types without consideration of the spatial variability in legacy P for given land use. The land use-based typical legacy P values were derived from the literature surveys [40,77] that compiled soil chemistry studies in the LOW. These databases indicated substantial variability in legacy P estimates for various land uses, likely stemming from variability in soils and historical changes in land cover and management. However, this variability remains unrepresented in WAM. For example, for improved pastures, the top horizon legacy P values ranged from 177 kg/ha to 730 kg/ha, while the calibrated TCNS model assigned a single value of 250 kg/ha to all improved pasture model cells. Similarly, legacy P values for active dairies ranged from 304 kg/ha to 9751 kg/ha, while a single calibrated value of 700 kg/ha was used. In order to illustrate the potential impact of this, the Baseline model accuracy was assessed at three upstream sub-basin locations in the TCNS (USGS02274010, USGS02274490, USGS02274505) (Figure 1), which were not included in the original model calibration-validation exercise. While the flow simulation performance of WAM at all three upstream locations was satisfactory, the same was not true for TP loads (Table S6). Simulated TP loads at USGS02274010 and USGS0224505 over-predicted observations by 41.6% and 179.4%, respectively, while at USGS02274490 TP loads deviated by 13.5% from observations. The worsened model performance in simulating TP loads at these upstream locations can be attributed to unaccounted spatial variability in legacy P values for the dominant land uses in these sub-basins (improved pastures, abandoned and active dairies) and is concerning as simulations showed that areas within these sub-basins had a very high UATPL (Figure 7).
The above analysis emphasizes the importance of accounting for spatial variability in legacy P while applying hydrologic and water quality models to identify critical sources, formulate and optimize nutrient control projects, and inform decision-making for watersheds with known legacy P problems. This is even more important when models are known to be sensitive to legacy P-related parameters. In order to achieve this, spatially accurate legacy P information is necessary. Recent advancements in chemical analysis techniques that accurately estimate soil phosphorus retention characteristics are an attractive alternative [95]. A comprehensive soil chemistry analysis project followed by geostatistical analysis techniques would be a viable approach in developing spatially accurate legacy P maps, as it would integrate effects of historical changes in land use and management as well as soil characteristics. Conducting model performance assessments at a finer scale (field, basin, or sub-basin) instead of watershed or a sub-watershed level would provide insights for optimizing a chemical analysis-based legacy P determination program and identifying any knowledge gaps in local P transport processes.

Implications for the Lake Okeechobee Watershed
Our results indicate that legacy P is responsible for 61-65% of the existing loads from the TCNS, implying its overall importance and the need for its specific consideration while identifying critical sub-regions, formulating solutions, and prioritizing projects; for example, the Targeted Restoration Area approach adopted in the Lake Okeechobee BMAP [61]. The same is true for other non-inorganic fertilizer P sources, which were found to contribute over 30% of 111.3 mt/y annual TP loads from the TCNS. One should be cautious in extrapolating relative source contributions to other LOW sub-watersheds. However, it is noteworthy that LOW sub-watersheds adjacent to the TCNS (Indian Prairie, Lower Kissimmee, and Fisheating Creek) also have agriculture as the dominant land use (over 50% in each), with 44% to 60% of their drainage areas being classified as pastures, the majority of which are improved pastures [101]. This suggests that the BMPs tailored to reduce legacy P and improve manure management/applications would be beneficial for those sub-watersheds as well. Nair et al. [102] recommended the adoption of a "soil phosphorus storage capacity" to determine the site-specific phosphorus fertilization rates to first exhaust legacy P via crop uptake without impacting productivity. This study in Florida showed promising results as a viable way of mining legacy P through reduced fertilizer applications but could benefit from further research and development. Pavinato et al. [99] discussed other practices from the perspective of legacy P, which includes soil amendments, crop rotation, efficient crop varieties, modern fertilizers, and inoculation with P solubilizing microbes, etc. Their suitability for the LOW conditions needs assessment, followed by incorporation into agricultural and urban BMPs programs. The results of this study highlight the need for expediting projects and programs to help further mine legacy P that exists on the landscape. Projects such as the water quality monitoring network enhancement, the innovative treatment technology application for P removal in the S191 basin [103], and the recently announced new STA [104] are steps in the right direction. Monitoring for the reduction of legacy P at the source level can be incorporated as a parameter to gauge the success of nutrient management BMPs. Even though the contributions of inorganic fertilizer management BMPs may seem low to overall TP load, these should not be interpreted as unimportant. These are crucial to prevent water quality from further deterioration and essential due to their low costs and potential of producing results in a relatively short timeframe when implemented properly.

Conclusions
This study presented a quantitative analysis estimating the contributions of various phosphorus sources (legacy P, inorganic fertilizers, and other P sources) to TP loads reaching Lake Okeechobee from the Taylor Creek Nubbin Slough sub-watershed using the Watershed Assessment Model. The relative importance of legacy P and fertilizer application rates was also assessed by performing a global sensitivity analysis. Legacy P was found to be the most dominant TP contributor and estimated to account for 63% of the existing loads from this sub-watershed. The sensitivity analysis results showed that under plausible changes in legacy P and fertilizer application rates, consistent with advanced fertilizer management, legacy P in improved pastures was the most influential model parameter affecting TP loads. The land uses that contributed the highest to TP loads were found to be improved pastures, active and abandoned dairies, citrus, and field crops. These results suggest that developing and implementing suitable conservation practices that can mine legacy P at the source level and expediting large-scale stormwater treatment area projects could meaningfully improve TP loads originating from legacy P from entering Lake Okeechobee. Future efforts should include the evaluation of impacts of climate change on legacy P dynamics in this region.
This study is a rare example of a mechanistic hydrologic and water quality model application estimating the impacts of legacy P on water quality. Our results are in agreement with other studies that demonstrated the importance of apportioning contributions of different phosphorus sources. A crucial need to develop spatially accurate legacy P information, and to incorporate that information into modeling frameworks for improved and reliable identification of critical source areas in watersheds that are likely to exhibit legacy P, was recognized in this work.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/land10090977/s1: Table S1: Summary of land use coverage in the TCNS basins (expressed as % area of each basin). Table S2: Summary of physical and chemical properties of soils within the TCNS sub-watershed. Table S3: Summary of Baseline WAM TCNS setup Performance During the Combined Calibration-Validation Period. Figure S1: Comparison of observed and WAM simulated outputs over the calibration-validation period in the TCNS basins for (a) S191 flow, (b) S191 TP load, (c) S154 flow, (d) S154 TP load, (e) S133 flow, (f) S133 TP load, (g) S135 flow, and (h) S135 TP load. Figure S2: Distribution of TCNS flows and TP loads among basins. Average annual values were from the 21-year (1993-2013) simulation period. Figure S3: Average annual TP load distributions under the Baseline and eight alternative scenarios at individual TCNS basins (a) S191, (b) S154, (c) S133, (d) S135, and (e) S154C. The Everglades Agricultural Area Model (EAAMOD) Overview. Table S4: Summary of legacy and inorganic fertilizer phosphorus parameters for the TCNS land uses under the Baseline and Alternative Scenarios. Table S5: Baseline WAM phosphorus attenuation parameters for land uses in the TCNS. Table S6