Nonlinear water clarity trends and impacts on littoral area in Minnesota lakes

Lake water clarity is an indicator of water quality, trophic status, and habitat condition. Changes in clarity impact lake ecosystems and may reflect land use changes or presence of invasive species. Quantifying temporal changes in water clarity can be challenging because clarity varies seasonally, annually, and spatially within and among lakes. We developed a hierarchical generalized additive model to quantify trends in water clarity (Secchi depth) from 1979 to 2018 for 909 Minnesota lakes, accounting for seasonal and spatial variability. Water clarity increased by 0.41 m across lakes from 1984 to 1988 and 2014 to 2018. Lake‐specific clarity trends varied: clarity did not change significantly in 59.0% of lakes, increased in 34.5% of lakes, and decreased in 6.5% of lakes. Water clarity dynamics caused considerable variability in littoral area between seasons and years. Our results have wide applications in aquatic ecology, including understanding changes to food webs, assessing fish habitat, and evaluating impacts of invasive species.

, and optical conditions (Hansen et al. 2019). Water clarity is also widely monitored as an index of trophic state and water quality in lakes (Carlson 1977;Farnaz Nojavan et al. 2019). Thus, the ability to assess changes in water clarity over time in individual lakes is crucial for understanding and predicting changes to populations and species interactions, managing lake organisms based on their habitat requirements, and ensuring that lakes are meeting water quality standards.
Water clarity impacts lake ecology and food webs through many pathways. For instance, clarity and morphometry determine the extent of a lake's littoral zone (Vadeboncoeur et al. 2008). The size of the littoral zone determines the relative benthic vs. pelagic contribution to whole lake primary production (Vadeboncoeur et al. 2002). Consumer diets switch from periphyton to phytoplankton as water clarity decreases and the littoral zone diminishes (Vadeboncoeur et al. 2003). Some organisms also utilize littoral resources more in some seasons than others (Hayden et al. 2014;Stewart et al. 2017). Understanding how seasonal and annual fluctuations in water clarity influence lake littoral area can shed light on food web changes and seasonal habitat requirements.
Lake water clarity may change over time for a variety of reasons. Clarity may decline in lakes due to algal blooms caused by increased nutrient inputs (Carpenter et al. 1999;Schindler 2006). Brownification due to increased terrestrially derived dissolved organic carbon (DOC) inputs reduces water clarity and is connected to land use change (Kritzberg 2017), increased precipitation (de Wit et al. 2016), and reduced acid rain (Evans et al. 2006). Increased precipitation can have variable effects on water clarity depending on trophic status, driving increased clarity in some eutrophic lakes due to dilution of phosphorus, decreased water residence time, and more stable stratification (Lisi and Hein 2019). Conversely, high water levels reduce the chance of winter kill of benthivorous fish and may lead to the loss of submerged macrophytes in shallow lakes, promoting a turbid-water state (Scheffer 2004). Water clarity may also respond to changes in biological communities, such as the increased clarity associated with invasion of zebra mussels (Dreissena polymorpha; Strayer 2009). In addition, lakes may become clearer over time following efforts to reduce external or internal nutrient loading (Jeppesen et al. 2007;Zamparas and Zacharias 2014).
Previous studies examining water clarity trends in Midwestern lakes have employed linear models to represent change over time, focusing on summer Secchi depth (Peckham and Lillesand 2006;Lottig et al. 2014;Olmanson et al. 2014). These studies suggest that Secchi depth has remained relatively stable in the majority of Midwestern lakes, and that individual lakes may exhibit significant water clarity trends in either direction. However, water clarity can be dynamic in response to the environmental and biological pressures outlined above, and linear models may not adequately capture these complex temporal trends. In addition, water clarity can exhibit substantial within-year temporal variability due to seasonal changes in light, temperature, and watershed inputs (Stadelmann et al. 2001). These seasonal trends are influenced by a lake's trophic status and timing of spring warming and fall cooling. A lake may also show considerable spatial variability in water clarity at different lake survey sites. Thus, summarizing Secchi depth by averaging observations across seasons or sites can muddy analysis of annual temporal trends, especially if the monitored dates or sites differ between years. These problems can be handled by restricting Secchi observations to one season (the approach utilized in the clarity trend studies above) or one site within a lake, but this approach results in lost information and effort when data are available for multiple seasons and sites. A model for assessing water clarity is needed that can account for nonlinear change within and between years, as well as spatial hierarchy within lakes.
We developed a hierarchical generalized additive model (HGAM) that estimates lake-and site-specific daily Secchi depth over multiple years. Our objectives were to quantify trends in water clarity for Minnesota lakes both within and among years, and to attribute heterogeneity in seasonal trends to differences in lake trophic status and timing of seasonal warming and cooling events. Our model improves upon previous approaches for analyzing trends in lake water clarity by estimating nonlinear seasonal and annual trends, predicting daily lake-specific estimates of Secchi depth, and using all available Secchi data for a lake regardless of season and site. We examined seasonal and annual clarity trends for 909 Minnesota lakes from 1979 to 2018. We also show that water clarity dynamics impose considerable variability in derived lake habitat metrics, using littoral area as an example. Our results have wide applications in lake ecology and management, including understanding seasonal and annual changes to food webs, assessing optical and thermal fish habitat, evaluating impacts of invasive species, and monitoring lake water quality.

Secchi data and covariates
Our study included lakes over 2 ha in Minnesota, USA. We compiled Secchi depth records through 2018 from the Minnesota Pollution Control Agency and the Minnesota Department of Natural Resources (MN DNR; Vitense and Hansen 2021). If data were available for sub-basins of multi-basin lakes, we kept the data at the sub-basin scale whenever possible. If > 1 Secchi observation occurred on the same date in a lake or sub-basin, we used the daily median. To account for different seasonal clarity trends in lakes with differing sources of turbidity, we compiled remotely sensed colored dissolved organic matter (CDOM) measurements for all available lakes (absorption coefficients at 440 nm, a 440 , averaged over years 2015-2016) from Olmanson et al. (2020). For analysis of long-term trends in water clarity, we selected lakes with Secchi records in ≥ 10 yr since 1979, ≥ 1 observation in each decade , and CDOM measurements. We also included Red Lake, for which observations began in 1990. This resulted in 909 lakes suitable for time series analysis (1008 time series total including subbasins) and 305,803 daily Secchi depth records. We used the CDOM data to classify these lakes/sub-basins as "algae-dominated" (a 440 < 3.38; n = 851) or "stained" (a 440 ≥ 3.38; n = 157; Supporting Information). We used the R package geoknife  to download gridded daily air temperature data associated with the spatial location of each lake from 1979 to 2018 (NARR data provided by NOAA/OAR/ ESRL PSL, Boulder, Colorado, USA). We used air temperature data to compute accumulated degree days above 0 C (DD0) at the date and location of each Secchi observation to account for annual differences in seasonal water clarity patterns attributable to temperature.

Secchi model
We developed an HGAM to model water clarity time series for all lakes and sub-basins. We assumed Secchi depth of subbasin i of lake l on day d of year y (S i,l,d,y ) followed a Gamma distribution, where basin of a lake that has > 1 sub-basins represented in the dataset (I Sub i ¼ 0 otherwise), C i {algae-dominated, stained} denotes whether the sub-basin was classified as "algae-dominated" or "stained," otherwise), and f ðÞ denotes a smoothed function of the predictor variable(s). The function f DOY i,d , DD0 i,d,y À Á is a tensor product smooth between day of year (DOY) and degree days (DD0); we used a cyclic cubic regression spline for the DOY effect and a thin plate regression spline (TPRS) for the DD0 effect. This function interacts with C i such that a different smoothed interaction is fit for algae-dominated vs. stained lakes/sub-basins. The function f Yr y À Á is a global year effect across all lakes, which we fitted using a TPRS. The function f Yr y ,Lk l À Á is a factor-smooth interaction between year and lake, where different year smoothers are fit for each lake, but all lakes share the same smoothing parameter (Pedersen et al. 2019). We used a TPRS for each lake-specific smoother. Equation 2 denotes that the variance of each Secchi observation is proportional to its expected value E S i,l,d,y À Á À Á . Supporting Information contains additional information and justification regarding model structure.
We fit the model using the bam() function from the mgcv R package (Wood et al. 2015;Pedersen et al. 2019;Li and Wood 2020). We employed the fast REML method and functionality to discretize covariates to speed up computation time. We used the function gam.check() to check model fit and ensure the basis dimension (k) for each smooth was adequate. We performed 10-fold cross validation using stratified partitioning to divide the data into 10 groups such that each lake's observations were divided roughly evenly among the partitions. We refit the model 10 times, leaving each partition out in turn, and used each fitted model to predict the Secchi depths of the data in the partition that was left out. We assessed predictive accuracy by computing the mean absolute error across all observations and seasonal errors for spring (April-May), summer (June-August), fall (September-October), and winter (November-March).
We used two complementary approaches to summarize temporal trends. To assess long-term change in water clarity across Minnesota, we computed differences in 5-yr averages for the latest years in the time series (2014-2018) and 30 yr prior (1984)(1985)(1986)(1987)(1988) for both the global trend (lake random effects set to zero, other predictors held at median values: DOY = 202, DD0 = 1743 C, algae-dominated lake) and lake-specific trends (lake-specific year smooths, intercepts, and algae-dominated/ stained class, with DD0 fixed at each lake's median for DOY = 202). We computed 95% confidence intervals for these differences using posterior simulation of the model parameters (Simpson 2018) and considered lake-specific change in Secchi depth to be statistically significant if the confidence interval did not contain 0. We also identified years of significant change in the time series for both the global trend and individual lakes by computing derivatives with respect to year using the method of finite differences and estimating 95% simultaneous confidence intervals for the derivatives at each year using posterior simulation (Simpson 2018). We defined years of significant change as years where the simultaneous confidence interval did not contain 0.

Seasonal and annual littoral area estimates
Using the fitted HGAM, we generated daily Secchi depth predictions for all lakes and sub-basins from 1979 to 2018. We converted these daily Secchi depth Z SD ð Þ estimates to daily light extinction coefficients K d ð Þ using the equation K d ¼ 1:7=Z SD (Poole and Atkins 1929). We then estimated the daily depth of the euphotic zone Z EZ ð Þ, defined as the depth at which photosynthetically active radiation (PAR) is 1% of radiation at the lake surface (Schindler 1971): Of the 1008 lakes and sub-basins with K d times series, 908 had hypsographic data such that we could calculate the area of the littoral zone, which we defined as lake surface area for which lake depth was equal to or shallower than Z EZ . We standardized lake littoral area by dividing littoral area by total lake area for comparison between lakes. We computed annual littoral proportion as the median daily littoral proportion in the ice-free seasons (spring-fall) of each lake-year, and seasonal littoral proportions were computed as the median daily littoral proportion within spring, summer, and fall of each year. We also computed median seasonal and annual littoral zone proportion across all lake-years for each lake.

Secchi model
The Secchi depth model passed diagnostic tests, and the deviance explained was 79.2% (note this is a 15.7% increase in deviance explained over a model that accounts for only differences in mean lake clarity). The overall cross-validated mean absolute error was 0.61 m; the seasonal mean absolute errors were 1.02 m for spring (n = 37,872), 0.57 m for summer (n = 206,709), 0.49 m for fall (n = 59,678), and 1.08 m for winter (n = 1541). Little bias was observed overall (0.02 m, where positive values indicate predictions were clearer on average than observed) and for summer and fall (0.00 and À0.05 m on average, respectively). Secchi depth predictions were biased slightly high for spring and winter (0.24 and 0.16 m on average, respectively). In addition, some lakes have high uncertainty in Secchi depth estimates at the tails of the time series where data are lacking (e.g., Fig. 1D).
A typical algae-dominated lake exhibited high water clarity in winter, reduced clarity during spring turnover, a spring clearwater phase, low clarity in summer and fall, and an increase in clarity when the season changed back to winter (Fig. 1A,B). Temperature interacted with day of year to influence water clarity; for example, lakes were more likely to have a clear-water phase as DD0 increased in the spring, and water clarity decreased as DD0 increased in the summer and fall. Stained lakes also exhibited higher clarity in spring/winter than summer/fall, but they had lower water clarity overall and less seasonal variation compared to algae-dominated lakes (Fig. 1C,D).

Water clarity trends
Water clarity increased by 0.41 m (95% CI: À0.15, 0.97) for a typical lake from 1984 to 1988 and 2014 to 2018 (i.e., when lake random effects are set to zero and other predictors are held at median values). However, lake-specific change between these periods at the ends of the time series varied across the state ( Fig. 2A). Average Secchi depth did not change substantially or was highly uncertain in 59.0% of lakes, while water clarity significantly increased in 34.5% and decreased in 6.5% of lakes (Fig. 2B). The global increase in clarity across all lakes occurred fairly gradually over the time series , with the longest span of stable water clarity occurring between 1994 and 2005 (Fig. 2C). Years with significant positive rates of change in the global trend generally correspond to years with a greater proportion of lakes with increasing trends than decreasing trends; similarly, more lakes had decreasing trends than increasing trends in years for which the global trend had significant negative rates of change. In addition, some lakes exhibited sharp increases or decreases in clarity over the course of the time series. For example, in Lake Geneva, water clarity increased substantially following invasion by zebra mussels in 2009, after a period of relatively stable clarity between 1979 and 2008 (Fig. 3A,B).

Littoral area
Lake littoral area varied substantially across lakes, years, and seasons in response to Secchi depth dynamics and lake morphometry. The median annual littoral zone proportion ranged from 0.06 to 1.00 across lakes from 1979 to 2018 (median littoral proportion = 0.75; Fig. 4A). Within lakes, littoral area was variable among years; within-lake inner 95% ranges (i.e., difference between 97.5 th and 2.5 th quantiles) of annual littoral proportion varied from 0.00 to 0.96 (median range = 0.23; Fig. 4A).
Lake littoral area varied between seasons due to seasonal fluctuations in water clarity. Littoral proportion was generally higher in the spring when water tends to be clearer: median spring littoral proportions were À0.03 lower to 0.58 higher than median annual littoral proportions (median difference = 0.05; Fig. 4B). Littoral area was reduced in summer and fall: summer medians were À0.41 lower to 0.02 higher than annual medians (median difference = À0.03), and fall medians were À0.40 lower to 0.00 higher than annual means (median difference = À0.04).
The above points are further demonstrated by examining how lake littoral proportion is influenced by changing Secchi depth for Mille Lacs, Minnesota (Fig. 4C). The littoral proportion of Mille Lacs increases dramatically as Secchi depth rises to 4.5 m, above which the entire lake is littoral. Since 1980, the littoral zone of Mille Lacs has generally expanded due to increased clarity. Mille Lacs has a larger littoral zone in the spring than summer and fall (median spring littoral proportion is 0.52; median summer and fall littoral proportions are 0.28 and 0.29, respectively). One lake with an estimated change of +13.5 m was excluded from the histogram. (C) Global long-term trend across all lakes (estimates shown for an algae-dominated lake, DOY = 202, DD0 = 1743 C, lake random effects set to zero) with 95% simultaneous confidence interval (dashed lines). Blue points denote years in which clarity was significantly increasing, and brown points denote years in which clarity was significantly decreasing based on a 95% simultaneous confidence interval for the derivative of the year smooth. (D) Proportion of lakes in each year that had nonsignificant rates of change, significant positive rates of change (increasing clarity), and significant negative rates of change (decreasing clarity) based on 95% simultaneous confidence intervals for each lake.

Discussion
We developed an HGAM for analyzing lake clarity trends, improving upon previous approaches by accounting for nonlinear seasonal and annual clarity dynamics. Approximately 1/3 of our study lakes had higher clarity in 2014-2018 compared to 1984-1988, only 1/15 had lower clarity, and 3/5 of lakes either had little change in clarity and/or the estimated change was highly uncertain. These trends are similar to those found for 122 Minnesota lakes from 1970 to 1992, where 53% exhibited no change, 36% significantly increased, and 11% significantly decreased in clarity (Heiskary et al. 1994). Our results are also similar to Lottig et al. (2014), who found that Secchi depth increased overall in Midwestern lakes from 1938 to 2012, with more lakes increasing than decreasing in clarity (7% vs. 4%, respectively). Our results differ somewhat from a study of satellite-derived Secchi depths finding slightly more Minnesota lakes had declining clarity trends than increasing trends from 1985 to 2005 (6.2% vs. 4.6%, respectively) (Olmanson et al. 2014). However, our approach for examining temporal trends differed from these studies in that we accounted for Secchi data collected outside of summer, and we used nonlinear models to assess trends over time both globally and for individual lakes. We showed that clarity for the typical lake in Minnesota has increased overall since 1979 with occasional years of decreasing clarity.
Lake water clarity dynamics caused littoral area estimates based on 1% PAR to vary seasonally and annually. Analyzing changes in littoral area over time based on time-varying water clarity estimates could improve our understanding of longterm and seasonal changes to lake food webs (McMeans et al. 2015) and could be used to predict the persistence of macrophyte-dominated stable states in shallow lakes (Scheffer 2004). In addition, this variability in littoral habitat has implications for fisheries management. For example, MN DNR Fisheries stocks walleye at a rate determined by lake littoral area, defined as surface area less than 15 ft (https://www. dnr.state.mn.us/fisheries/management/stock.html). However, littoral estimates based on clarity are often substantially confidence intervals for daily mean estimates. (B) Long-term trend for Lake Geneva (at DOY = 202 and DD0 = 1766 C, the lake's median DD0 at DOY = 202) with 95% simultaneous confidence interval (dashed lines). Blue points denote years in which clarity was significantly increasing, and brown points denote years in which clarity was significantly decreasing based on a 95% simultaneous confidence interval for the derivative of the lake's year smooth. Lake Geneva was listed as infested with zebra mussels in 2009. different than those based on depth alone (Fig. 4C), where estimates based on depth tend to estimate less littoral area for lakes with high clarity and more littoral area for lakes with low clarity compared to estimates based on Secchi depth (Fig. 5).
Our approach is useful for evaluating seasonal and annual clarity patterns for many lakes. Importantly, although many of the lakes in this dataset do not have year-round Secchi depth observations to enable estimation of lake-specific seasonal trends, our landscape-scale model allows information about seasonal patterns to be shared among lakes based on lake color. The algae-dominated pattern was similar to the broad Midwestern pattern observed by Read et al. (2017) and is representative of the majority of lakes in our study. Stained lakes had more muted seasonal dynamics, consistent with lower algal growth and total suspended solids in lakes with high CDOM (Brezonik et al. 2019). In addition, the model hierarchy allows information to be shared across sub-basins in multibasin lakes to estimate a lake-wide temporal trend, while accounting for differences in clarity between sub-basins. Our model also accounts for season when looking at long-term changes in clarity over time, enabling use of all data regardless of season. We emphasize again that estimation of nonlinear temporal trends allows for a more nuanced investigation of change over time, including identification of years of  significant change globally and for individual lakes. Lakespecific temporal trends are useful for evaluating the impacts of in-lake management actions, invasion of zebra mussels, or land use changes. Finally, modeled Secchi depth time series can be used to generate time series and associated uncertainty for derived habitat metrics, such as optical and thermal fish habitat, with broad applications in lake ecology and management.