Windows into the past: lake sediment phosphorus trajectories act as integrated archives of watershed disturbance legacies over centennial scales

Historic land alterations and agricultural intensification have resulted in legacy phosphorus (P) accumulations within lakes and reservoirs. Internal loading from such legacy stores can be a major driver of future water quality degradation. Yet, little is known about the magnitude and spatial patterns of legacy P accumulation in lentic systems, and how watershed disturbance trajectories drive these patterns. Here, we used a meta-analysis of 113 paleolimnological studies across 124 lakes and four reservoirs (referred here on as lakes) in 20 countries to quantify the linkages between the 100 year trajectories of P concentrations in lake sediments, watershed inputs, and lake morphology. We find five distinct clusters for lake sediment P trajectories, with lakes in the developing and developed world showing distinctly different patterns. Lakes in the developed world (Europe and North America) with early agricultural intensification had the highest sediment P concentrations (1176–1628 mg kg−1), with a peak between the 1970–1980s and a decline since then, while lakes in the developing world, specifically China, documented monotonically increasing sediment P concentrations (857–1603 mg kg−1). Sediment P trajectories reflected watershed disturbance patterns and were driven by a combination of anthropogenic drivers (fertilizer input and population density) and lake morphology (watershed to lake area ratio). Specifically, we found the largest legacy accumulation rates to occur in shallow lakes experiencing long-term land-use disturbances. These links between land-use change and P accumulation in lentic systems can provide insights about inland water quality response and help to develop robust predictive models useful for resource managers and decision-makers.


Introduction
Exponential increase in population and food demands in the past century have led to agricultural intensification and land use change (Tilman et al 2011). Long term changes in cropping patterns, excessive fertilizer applications, and livestock operations have shifted the global phosphorus (P) cycle outside its natural historic range of variability (Vollenweider 1971) and contributed to eutrophication and hypoxia in inland and coastal waters. Algal blooms and hypoxia is a long standing concern in lakes around the world (Schindler 2012, Michalak et al 2013, Nürnberg et al 2013, Ho et al 2019.
Decades of efforts to manage watershed P loading have yielded limited results. Despite reduction in total P loads observed in few systems (Sharpley et al 1996, Baker et al 2014, most lake P concentrations seem to have either attained stasis or continued to increase in various regions around the world , Oliver et al 2017. It is now recognized that this lack of success in water quality improvement can be partially attributed to the build-up of excess P in soils, groundwater, rivers, and lake sediments over decades of agricultural intensification and increasing population densities (MacDonald and Bennett 2009, Keatley et al 2011, Haygarth et al 2014, Goyette et al 2016, Powers et al 2016, Stackpoole et al 2019, van Meter et al 2021. These legacies contribute to elevated concentrations and often decadal lags in ecosystem improvements even when inputs have ceased , van Meter et al 2018, Basu et al 2022. Lake sedimentary deposits can serve as integrated archives of P accumulation, and watershed disturbance, and in turn provide insights on the deviations from baseline ecosystem conditions and the temporal trends in water quality degradation (Smol 1992, Dearing et al 2008, Bhattacharya et al 2016, Edlund et al 2017. Such long-term perspectives can be a critical benchmark for validating watershed nutrient transport and retention models (Anderson et al 2006, van Meter et al 2017, and for effective management of inland and coastal aquatic ecosystems (van Meter et al 2018). Typically, studies of lake sediment P trajectories have focused on reconstructing ecological evolution of lakes using sediment stratigraphic records of ecological indicators (e.g. diatoms and algal pigments), along with radio-active isotopes (e.g. 210 Pb and 137 Cs) used to establish the chronological history (Appleby 2008). These studies have mostly focused on understanding the drivers of sediment P dynamics in individual lakes (Engstrom et al 2006, Otu et al 2011, Levine et al 2018. Regionalscale studies that explore P trajectories are limited, albeit a few exist, such as a study across two dozen lakes in Finland (Tammeorg et al 2018), or a range of shallow lakes in Southern Florida, USA (Kenney et al 2016). Using a meta-analysis approach, Carey and Rydin (2011) show that P burial patterns are representative of lake trophic status, alluding to the linkage between lake sediment and watershed inputs. Keatley et al (2011) used trajectories of sediment diatom assemblages to reconstruct lake water P concentrations and relate them to watershed disturbances. The study largely focused on the lakes in North America and Europe, and further warrants an exploration of sediment P trajectories across a wider geographical area.
Here, we build on this body of work and use a meta-analysis approach to explore century scale trajectories in lake sediment P across the world, and relate these trajectories to watershed disturbance legacies. First, we use a machine learning approach to identify dominant clusters of sediment P trajectories, and analyze the relationship between the temporal trajectories of the clusters and the corresponding land use trajectories at the cluster level. The cluster level analysis is complemented by lake-level analysis where we explore the relationships between the total P accumulation in the sediment, and cumulative P inputs from the watershed. The cluster level analysis allows us to explore the temporal trajectories, while lake-level analysis allows us to better elucidate dominant controls on total accumulation.

Data synthesis for meta analysis
We used the following keywords 'paleolimnology' , 'phosphorus' , 'lake' , 'sediment' , 'core' to search the Scopus database (2020) for publications related to lake sediment cores. We applied additional filters based on subject areas (i.e. Earth and Planetary Science, Environmental Science, Biological and Agricultural Science and their related subject areas), document type (article), and language (English). Availability of 210 Pb chronology was used as additional key criteria for paper selection. Chronologically constrained sedimentary records were essential for our study as it allowed the comparison of sediment P trajectories with watershed disturbance legacies. Our Scopus database search resulted in a total of 600 papers, and additional screening for 210 Pb chronology led to final 200 papers. Then, we used the 'digitize' R package (Poisot et al 2016) to extract the sediment P trajectories and age-depth chronology from these papers to be used for further analyses.
The P trajectories for each sediment core was obtained by plotting the digitized depth-specific sediment P concentrations (P sed ; mgP kg −1 dry sediment) against its depth-chronology extracted from the publications. Age chronologies are developed using radiometric dating techniques (e.g. 210 Pb isotopes; Appleby 2008) that inform about the age/year at which the sediment P was deposited at the lake bottom. For lakes with multiple sediment cores, a basinwide single sediment P trajectory was estimated by aggregating the individual cores chronologically.

Lake and watershed attributes
Lake coordinates and mean lake depth were extracted from the published journal articles. Lake area was estimated from the lake polygon database Hydro-LAKES (Lehner et al 2008). We delineated the drainage areas of the lakes using the digital elevation model from the global hydrographic database HydroSHEDS (Lehner et al 2008). We used six metrics to describe watershed disturbance: (a) C per (%) = crop percent, (b) P fert (kg ha −1 watershed area/year) = P fertilizer input per unit watershed area, (c) P fert n (kg ha yr −1 ) = P fert × W:L area , (d) N p (persons ha −1 ) = population density, and (e) N p n (persons ha −1 ) = N p × W:L area . The normalization of the fertilizer inputs and the population density by the watershed to lake area ratio allows us to integrate both watershed inputs and lake morphology. The annual trajectory  of crop percentage (C per ) was obtained from a historical cropland dataset (Ramankutty and Foley 1999) aggregated to the watershed scale by multiplying the percent crop to watershed area. P fertilizer trajectories  were estimated from the dataset of P fertilizer use on croplands (kg ha −1 cropland yr −1 ) (Lu and Tian 2016), and the corresponding cropped area of each watershed. In addition to fertilizers, manure application to croplands can be an important source of P to the lakes, though fertilizer still constitutes between 65% and 70% of the inputs (Bouwman et al 2017, van Puijenbroek et al 2019, Glibert 2020). Current long term data availability limits the use of the manure dataset, but this should be explored in future work. Finally, gridded population density trajectories  were obtained from the HYDE 3.2 dataset (Goldewijk et al 2010). It is important to note that these reconstructed data products, though widely used, have uncertainties that might confound our interpretations; however, these are some of the most widely used datasets and the uncertainties are much less in the more recent timeframe (1960s to current).

Statistical methods
To explore the spatio-temporal variability in lake sediment P trajectories, we utilized a two-pronged approach. First, we used a hierarchical clustering method to identify the major groups of lakes with similar temporal patterns (section 2.3.1). The median sediment P concentration trajectories for each cluster were then compared with the potential drivers described above in section 2.2. The cluster level analysis was complemented by a lake-level analysis where we analyzed the dominant controls of lake-sediment P accumulation using ordination and regression methods (section 2.3.2). The lake level analysis presents temporally aggregated metrics and allows us to explore spatial variability and dominant controls. In contrast, the cluster level analysis is aggregated spatially but allows us to explore temporal trajectories.

Hierarchical clustering of sediment P trajectories and changepoint analyses
Agglomerative hierarchical clustering-an unsupervised machine learning algorithm that utilizes a 'bottoms up approach' (Kaufman and Rousseeuw 2005), was used to find lake groups with similar lake sediment P trajectories across all our study lakes (n = 128). Hierarchical clustering is a well established, and robust technique that does not require predetermined number of clusters as an initial parameter, and can work with time series data with varying levels of similarities, which makes the method more appropriate for the real world time-series data (Aghabozorgi et al 2015). As such, hierarchical clustering has been successfully used for clustering long-term water quality and sediment time series (e.g. (Astel et al 2006, Legendre et al 2012, Dugan et al 2017, Byrnes et al 2020. To allow for cross-site comparison, we developed z-score trajectories by standardizing P sed values in each core to a mean of zero and variance of one. The z-score trajectories were smoothed using a moving average (span = 3) (TTR package v. 0.24.2, Ulrich 2016) to remove data noise or abrupt inter-annual variability (Bigler and Hall 2003). The smoothed z-score trajectories were clustered using Ward's linkage hierarchical clustering method based on squared Euclidean distances that minimises within group dissimilarity at each linkage (Ward 1963). We explored P sed trajectories between 1900 (±5 years) to 2010 (±5 years), which resulted in some trajectories of unequal length. To account for unequal trajectories and an additional check for our clusters, we used the dynamic time warping as the distance measure in our hierarchical clustering (Aghabozorgi et al 2015). All the clustering outcomes were evaluated by observing the cluster dendrograms and using the sum of squared errors (SSE) between clusters to test the cluster 'coherence' , with better clusters having low SSE (Aghabozorgi et al 2015). Finally, changepoint analysis was also conducted for each trajectory (changepoint package v. 2.2.2, Killick and Eckley 2014), to statistically detect any prominent shift in P sed trajectories.

Exploring lake level variability in sediment P trajectories and their controls
To understand spatial controls on sediment P trajectories at an individual lake level for all the study lakes (n = 128), we used the Principal Component Analysis (PCA), an ordination based analysis technique (Ter Braak and Šmilauer 2002). We used two metrics to integrate across the temporal trajectory in sediment P: (a) the rate of sediment P accumulation (P R ; mg kg yr −1 ) calculated as the difference in P concentrations between 1900 (±5 years)-2007 (±5 years) divided by the duration of P trajectory, and the (b) cumulative sediment P concentration totals over the duration of 1960 (±5 years)-2007 (±5 years) (P M ; mg kg −1 ). These two metrics were used to assess the spatial variability in P trajectories and allow comparison with watershed and lake variables, including lake surface area (L area ; km 2 ), ratio of watershed to lake area (W:L area ), average lake depth (D; m), and watershed disturbance indicators averaged over the duration of 1960-2007, including crop percent (C per ; %), and population density (N p n ; persons ha −1 ). We also used the cumulative fertilizer inputs (P fert ; kg) and cumulative normalized fertilizer input (P fert n ), as the summation of the P fert and P fert n metrics (section 3.4) over the timeframe 1960-2007.

Metadata synthesis
Our final database consisted of 124 lakes and four reservoirs (from here on both lakes and reservoirs will be referred to as lakes) spanning across 20 countries and seven continents (figure 1(A); table S1 available online at stacks.iop.org/ERL/17/034005/mmedia). Largest number of lakes in our study were in China (n = 43), followed by the USA (n = 37), UK (n = 23), Canada (n = 6), and Japan (n = 5). The rest of the 14 countries had less than 5 lakes in each. As such, the lakes located between 20 ′ N and 60 ′ N are over-represented, while the arctic and tropics are under-represented in our dataset. The distribution of lake surface area (L area ; km 2 ), watershed-lake area ratio (W:L area ), lake depth (D; m), and %watershed that is cropped (C per ; %) (figures 1(B)-(E)) across our study sites represent a gradient of lake morphological and watershed characteristics. In terms of lake morphology, 50% of the sampled lakes were shallow (<5 m), and another 25% of the lakes had depths ranging between 5 and 10 m, while 25% were deeper lakes ( figure 1(D)). The surface area spanned between 0.022-6.8 × 10 4 km 2 , with 40% of the lakes having a surface area <10 km 2 ( figure 1(B)). Land use across the lake basins varied significantly (C per ranged from 0.3% to 85%), with 50% of lakes located in pristine watersheds (C per < 20%), 15% of lakes in agriculture-dominated watersheds (C per > 50%), while in the remaining 35% of the lakes cropped area varied between 20% and 50% (figure 1(E)). Thus, our lake database over-represents pristine lakes in the northern hemisphere, with less number of studies from reservoirs and agriculturally impacted lakes in the developing world, and thus may create gaps in our understanding of global P accumulation.

Sediment P trajectories: cluster analysis
Sediment P concentration trajectories (P sed , mg P kg −1 ), spanning over the last 100 years, exhibit significant variability in space and time (figure 2). Five distinct clusters were identified that described the behavior of the 128 sediment P cores (figures 2(A)-(F); table S1). Cluster 1 (C1) includes 51 lakes, and is characterized by minimal temporal trend and no significant changepoint over the 100 year timeframe (figures 2(A) and (B)), and the lowest P sed ( figure 2(A)). Majority of the C1 lakes are located in the USA (51%), followed by China (37%), and Canada (7%) (figures 2(G) and 3). In contrast to the C1 cluster, the other four clusters all show a clearly defined changepoint after which there is a significant increase in P sed . The median changepoint for cores in cluster 2 (C2) is 1978, while it is 1946 for cluster 3 (C3), 1937 for cluster 4 (C4) and 1960 for cluster 5 (C5) (figures 2(B)-(F)). C2 lakes have a much later changepoint than the other clusters, and exhibit a monotonically increasing concentration trajectory beyond the changepoint. Lakes in this cluster are located mostly in China (68%), followed by the UK (18%), and USA (10%) (figures 2(G) and 3). Agricultural intensification occurred much later in China compared to the other regions of the world, and this possibly explains the later changepoint.
C3 and C4 lakes show a peak around the 1990s (figure 2(D)) and 1970s (figure 2(E)), respectively, and a declining trend after that. The peak in C4 lakes is earlier, while the decline is over a longer timeframe and more dramatic than in C3 lakes. Lakes in C3 (n = 26; figures 2(C), (D) and (G)) and C4 (n = 12; figures 2(E) and (G)) belong primarily to the developed world (C3: ∼40% lakes in USA and UK, 14% in China; C4: 63% in UK, in USA, and 13% in Japan) (figures 2(G) and 3). Finally, C5 included only 11 lakes spread mostly across the US (43%) and UK (29%), with few lakes in Finland. These lakes recorded a steep increase in sediment P concentrations from the 1960s till 2000s (figures 2(F) and 3).
We further find that the P sed between 1960 and 2010 in C3 (median: 19 210 mg kg −1 ) and C4 (median: 33 780 mg kg −1 ) lakes to be higher than in C2 (median: 12 340 mg kg −1 ) lakes, alluding to the higher degree of historical impact on lakes in the developed world. The high P sed and the distinct peak in the 1970s-1990s in C3 and C4 lakes are indicative of highly impacted systems that might be slowly recovering. Indeed, the peak and decline response in C3 and C4 lakes, is parallel to peak and decline response that was observed in the P surplus trajectories for the watersheds in the UK and USA (Powers et al 2016). C5 lakes also belong to the developed world, but in contrast to C3 and C4 lakes, P sed is monotonically increasing, indicating no distinct recovery behavior ( figure 2(A)).

Land use and P input trajectories for the clusters
To better understand the drivers of sediment P accumulation, we then compared the P sed trajectories (figure 2) with land use trajectories of their watersheds (figure 4). We found C1 lakes to have the lowest percent cropland (C per ), which comprised the majority of pristine lakes in our dataset, with no substantial change over time ( figure 4(A)), and this is consistent with the corresponding lack of changepoint in the P sed trajectories. There is also no significant trend in C per in both C3 and C4 lakes, although C4 lakes do show a decrease in percent cropland since the 1940s that correspond to the changepoint of 1937 for their P sed trajectories (figure 2(E)). This is in contrast to C2 and C5 lakes where C per increases over time, although the increases are of a much smaller magnitude than the corresponding P sed trajectories (figures 2(C) and (F)). It is further interesting to note that C2 lakes have the highest C per but lower P sed than C3 and C4 lakes (figures 2(A) and 4(A)). This suggests that land use trajectories are possibly not the best predictors of lake sediment P accumulation. Indeed, fertilizer application and population dynamics might be a better indicator of P sed , as they more accurately capture the lake inputs. For example, agricultural activities in China (C2 lakes) started early contributing to the high C per , but the extensive use of P fertilizers is a more recent phenomenon than in North America and Europe (C3 and C4 lakes), and this possibly is responsible for the more recent changepoint of the 1970s in the C2 lakes.
To explore the possible effects of watershed disturbance further, we compared P sed trajectories with P fertilizer use (P fert ; kg ha −1 watershed area/year) and population density (N P ; persons ha −1 ) trajectories in watersheds draining into the lake (figures 4(B) and (E)). We find C1 lakes to have the lowest P fert and N P ; this is consistent with their P sed trajectories ( figure 2(A)). In contrast, C2 lakes have a monotonically increasing P fert trajectory since the 1960s, while C3, C4 and C5 lakes have had a peak P fert between 1970 and 1990, and a declining trend since then. This corresponds well with the monotonically increasing P sed for the C2 lakes, and a peak and declining trend for the C3 and C4 lakes. Parallel also to the P sed trajectories the peak in P fert in C4 lakes precedes the peak in C3 lakes. For C5 lakes, even though P fert decreased after the 1980s, the continued increase in population density from 1960s-late 2000s likely contributed to the increase in sediment accumulation (figure 4).
While temporal patterns in P fert trajectories are similar to P sed trajectories the absolute magnitudes of P sed trajectories tell a different story. For example, while C2 lakes have the largest P fert , and N p , P sed in these lakes is significantly lower than C4 and C5 lakes (figure 2). Indeed, lakes with the highest P sed trajectories are C4, and C5 lakes followed by C3, C2 and C1 lakes, while the P fert magnitudes are the largest for the C2 lakes followed by C4, C3, C5 and C1 lakes. We hypothesize that this mismatch occurs because P fert ( figure 4(B)), which is the fertilizer inputs normalized to the watershed area, might not appropriately capture the interaction between the lake and its watershed. For example, the same magnitude of fertilizer inputs might impact a smaller lake differently than a larger lake.
To explore this further, we compared the ranges in the watershed to lake area ratio (W:L area ) between the various clusters. We found that C4 lakes have the highest W:L area (median: 71; IQR: 41-135), while C1 and C2 lakes have the smallest W:L area (for C2 median: 14; IQR: 8-43; figure S1). A higher W:L area indicates a larger watershed that contributes to greater inputs, and a much smaller lake area to buffer the larger input; this is possibly what leads to the higher P sed magnitudes. In fact previous studies on watershed and lake linkages have also established the unique impact of W:L area on nutrient loading and their subsequent retention within lakes (Schindler 1971). To explore the morphometric mediation of watershed disturbances on P sed , we compared the W:L area normalized fertilizer inputs and population (P fert n and N p n ) to the P sed trajectories (figures 4(C) and (F)). Consistent with the P sed trajectories, C4 lakes have higher magnitudes of P fert n compared to C2 lakes. C5 lakes are an exception with P sed among the highest and P fert n much lower; however it is important to note that C5 lakes have an extremely large range of W:L area ratios and the smallest number of lakes in its cluster, making it difficult to make cluster-specific generalizations (figure S1). To summarize, we find that lake P sed trajectories are driven by watershed inputs mediated by watershed and lake morphology, and thus the P fert n metric that explicitly considers both these controls is best able to describe the P sed trajectories.

Lake level assessment of sediment P magnitudes and accumulation rates
We then complemented the cluster-level analysis with lake-level analysis by exploring spatial patterns in sediment P accumulation magnitudes and rates (P M ; mg kg −1 and P R ; mg kg yr −1 ) across the world. Our samples were dominated by lakes in North America, Europe and China, with high accumulation magnitudes apparent in eastern US, UK and eastern China lakes ( figure 5(A)). The UK lakes have the highest P R (figure 5(B)) and P M (figure 5(C)), despite having lower fertilizer inputs (P fert ; cumulative fertilizer inputs from 1960 to 2007) than lakes in North America or China ( figure 5(D)). These are also the shallowest lakes, with the smallest depth ( figure 5(E)), and the largest watershed to lake area ratio (figure 5(F)), leading to the largest normalized fertilizer inputs ( figure 5(G)). These results are consistent with previous studies that report the link between external P loading and morphology on P accumulation in lakes. For instance, elevated accumulation of P sed in smaller lakes was also reported in Finnish (Tammeorg et al 2018) and other European lakes (Keatley et al 2011). Lakes with higher W:L area experience greater watershed inputs per unit lake area (Schindler 1971) that potentially leads to higher sediment P accumulation.
We further use the PCA to understand the drivers of P accumulation across all lakes ( figure 6(A)). We find that PCA axis 1 and 2 can together explain 53% of the variance of the P M and P R metrics. We found lake depth (D) and cumulative normalized fertilizer inputs (P fert n , cumulative normalized fertilizer inputs from 1960 to 2007) to emerge as two of the strongest controls on accumulation (figure 6(A)). The closer correspondence between P M andP fert n compared to P M andP fert is consistent with the cluster level analysis that argued for considering both watershed inputs and lake morphometry to explain lake sediment P accumulation (figure 6(A)). Indeed, while there exists a significant relationship between P M and P fert n (figure 6(C) and table S2), with higher P M in lakes with higher cumulative normalized fertilizer inputs (figure 6(B)), no significant relationship exists between P M andP fert (table S2). A significant negative relationship also exists between P M and D (table S2) with deeper lakes having smaller P M values (figure 6(D)). Depth emerges as a stronger control for the shallow UK lakes, with a long history of development, while cumulative fertilizer inputs, average cropped area, and population appear to have a greater influence for the lakes in China (figure 6(A)). Recent fertilizer intensive agricultural practices have been documented in several Asian countries, including China (Lu and Tian 2016) that have caused the sediment P enrichment (Rose et al 2011, Wu et al 2013. Our results highlight the differences in lake P M between the developed countries and developing countries, with developing country P M being more strongly driven by cumulative fertilizer inputs, while P M in developed countries (like UK) are driven by the morphometry (W:L area ) (figure 6(A)).
The lake level analysis thus confirms our cluster level analysis that found the normalized fertilizer application rate (P fert n ) to be the best predictor of sediment P accumulation. This metric takes into account both the cumulative fertilizer inputs to the lake (P fert ), and the lake morphometry that modulates the response, with smaller, shallower having greater productivity and thus accumulation magnitudes. The finding is important and is interesting and important and suggests that in absence of lake sediment core data, it might be possible to estimate lake sediment P accumulation magnitudes as a function of watershed P inputs and lake morphometry.  table S1 for details about all the study lakes.

Conclusion and implications
As inland waters are experiencing increased eutrophication and algal blooms, it is increasingly important to understand the role of lakes in mediating watershed scale nutrient export. Here, we use a meta analysis approach to synthesize lake sediment P concentration trajectories in 128 lakes across the world, and relate them watershed disturbance . The lakes were color coded based on their location. Countries with more than five study lakes (n > 5) in our database, namely the UK, USA, Canada, China, and Japan are indicated in the legend. Rest of the countries (n < 5 lakes) were shown by gray filled circles and the name of the country was added next to the symbol (B)-(D). Plots show the variability in average PM across ranges of cumulativeP fert , cumulative normalized P fert n , and depth (D).
legacies, specifically histories of agricultural intensification and urbanization over the last century. Using a machine learning approach, we identify five dominant lake sediment trajectory types, with distinct trajectories observed for the developed and developing regions of the world. We found 50% of the lakes sampled to be relatively pristine, with the lowest sediment P concentrations and no significant temporal trends (cluster 1 lakes). A second category of lakes had sediment P trajectories characterized by a peak in the 1970s to 1990s followed by a more recent decline (clusters 3 and 4). These lakes were characterized by the highest P concentrations, and are mostly located in the developed world (North America and Europe). In contrast, lakes in the developing world (China) had lower P concentrations in lake sediments, but an exponentially increasing accumulation pattern (cluster two lakes). The cluster level analysis was complemented by lake level PCA analysis to identify dominant drivers of lake sediment P accumulation. Both lake and cluster level analysis revealed the importance of watershed inputs (cumulative fertilizer inputs) and lake morphometry (depth and watershed to lake area ratio) on sediment P accumulation. Specifically, we observed a positive correlation between lake sediment P accumulation and cumulative fertilizer inputs to the watershed normalized by the watershed to lake area ratio, highlighting the combined role of watershed inputs and lake morphology on sediment accumulation patterns.
Our findings highlight the effect of wastewater management and agricultural development on the timing of sediment P accumulation and eutrophication across developing and developed countries. Fertilizer-driven increase in P sed in developed countries occurred earlier (1960s-1980s), and the current declining trends are driven by a decrease in cropped area and fertilizer applications, as well as detergent P bans and improvement of wastewater treatment practices that occurred across Europe and North America in the 1970s. These observations are typical for the developed world that had an early onset of watershed anthropogenic alterations followed by an attempt to manage and mitigate water quality issues in the 1970s-1980s (Müller et al 2007, Blumentritt et al 2013. In comparison, developing or in-transition countries, such as China, experienced a later onset of agricultural intensification, with extensive land reclamation for rice production, and the associated peaks in fertilizer application (Yousaf et al 2017). Despite the delay in agricultural intensification, the exponential rates of increase have led to rapid build-up of sediment P stores and have turned these regions into hotspots of eutrophication and coastal hypoxia (Bouwman et al 2009, Fink et al 2018.
Our study has three novel contributions that have implications for lake monitoring and management. First, our datasets highlight a focus on pristine lakes, and relatively fewer studies on reservoirs and more impacted lakes, which limits our ability to manage impacted systems. Further, the majority of study lakes are located in the northern hemisphere and developed countries, which creates gaps in our understanding of global P accumulation from the more agriculturally intensive, developing regions of the world. Second, our analysis highlights that lake sediments are effective integrators of watershed disturbance patterns, and can potentially act as a long term archive of watershed P loads. Indeed, given the lack of monitoring data over long time scales, it might be possible to utilize lake sediment data to calibrate coupled watershedlake models. Finally, we find that the spatio-temporal variability in lake sediment P accumulation is driven by a combination of anthropogenic drivers (fertilizer input and population density) and lake morphology (lake depth, watershed to lake area ratio), with shallow lakes in regions of intensive agriculture having the largest accumulation magnitudes. This highlights the role of past watershed P inputs that can potentially build up as legacy P in lakes, and act as a long term source of P through internal loading, fuelling eutrophication and algal blooms for years after watershed inputs have ceased (Jeppesen et al 2005, Orihel et al 2017, van Meter et al 2021. Such insights about inter-region variability in legacy P accumulation in lakes is critical for the development of robust predictive models that can be used for resource managers and decision makers for the management of aquatic resources.

Data availability statement
Papers used in the meta analysis are listed in the supplementary table S1.
The source for fertilizer usage, percent crops, and population data are included in the manuscript.
The data that support the findings of this study are available upon reasonable request from the authors.