Shifting Patterns of Summer Lake Color Phenology in Over 26,000 US Lakes.

Lakes are often defined by seasonal cycles. The seasonal timing, or phenology, of many lake processes are changing in response to human activities. However, long-term records exist for few lakes, and extrapolating patterns observed in these lakes to entire landscapes is exceedingly difficult using the limited number of available in situ observations. Limited landscape-level observations mean we do not know how common shifts in lake phenology are at macroscales. Here, we use a new remote sensing data set, LimnoSat-US, to analyze U.S. summer lake color phenology between 1984 and 2020 across more than 26,000 lakes. Our results show that summer lake color seasonality can be generalized into five distinct phenology groups that follow well-known patterns of phytoplankton succession. The frequency with which lakes transition from one phenology group to another is tied to lake and landscape level characteristics. Lakes with high inflows and low variation in their seasonal surface area are generally more stable, while lakes in areas with high interannual variations in climate and catchment population density show less stability. Our results reveal previously unexamined spatiotemporal patterns in lake seasonality and demonstrate the utility of LimnoSat-US, which, with over 22 million remote sensing observations of lakes, creates novel opportunities to examine changing lake ecosystems at a national scale.

Introduction

Lakes are critical freshwater resources that are highly sensitive to stressors such as climate change (Woolway et al., 2020) and altered land use (Martinuzzi et al., 2014). Globally, these stressors are shortening the duration of ice cover (Sharma et al., 2019), increasing rates of lake carbon burial (Heathcote & Downing, 2012), increasing evaporative water loss (Wang et al., 2018), warming surface waters (O’Reilly et al., 2015), and changing mixing regimes (Maberly et al., 2020), all of which influence lake productivity and ecological state. These changes manifest themselves in the seasonality of lake processes. Just like a deciduous forest that comes to life in the spring, inland water bodies are characterized by a predictable seasonal succession of biological processes (Sommer et al., 2012). In the spring, many lakes experience a diatom bloom, followed by a “clear‐water” phase where zooplankton rapidly devour the newly plentiful phytoplankton (Matsuzaki et al., 2020). Summer algal biomass is constrained by nutrient availability, with nutrient‐rich eutrophic lakes experiencing near‐constant summer phytoplankton blooms, and nutrient‐poor oligotrophic lakes experiencing relatively clear waters (Sommer et al., 1986). The difference between these states is visible to the naked eye, as the predominant color of a lake lies along a spectrum of blue (oligotrophic) to green (eutrophic); or as dissolved carbon concentrations increase, brown (dystrophic) (Webster et al., 2008).

The color of a lake reveals a lot about lake productivity and ecological state. A green lake will have a greater abundance of phytoplankton and a higher rate of carbon burial than a blue lake (Heathcote & Downing, 2012). Browning or greening of oligotrophic lakes may result in oxygen depletion and anoxic conditions (Knoll et al., 2018; Müller et al., 2012), which impacts nutrient cycling. Shifts in the magnitude and timing of annual color changes are indicators of short‐term external (weather, nutrient, and carbon loading) and internal (biology) factors and/or long‐term climate, watershed, and food web changes. These changes are not confined to single lakes, with landscape‐level drivers impacting the color regimes of entire regions. For instance, shortened ice cover durations (Sharma et al., 2019) are shifting the spring‐phytoplankton bloom earlier (Winder & Schindler, 2004), increases in dissolved organic carbon are browning lakes (Monteith et al., 2007; Roulet & Moore, 2006), and invasive zebra mussels are increasing water clarity (Binding et al., 2007), all at regional scales.

For a single lake, observing the annual pattern of lake color provides insight into the local ecosystem. At larger scales, simultaneously observing the annual patterns of many lakes provides evidence of the impacts of climate and land‐use change and is critical in understanding the role of inland waters in carbon production and sequestration. Remote sensing enables macroscale analysis because it captures a wide range of hydrologic conditions (e.g., Allen et al., 2020) with regular sampling intervals and global coverage. The Landsat series of satellites specifically provides over three decades of observations and can be used to accurately estimate water quality parameters such as chlorophyll‐a (Cao et al., 2020; Lin et al., 2018), colored dissolved organic matter (CDOM) (Griffin et al., 2018; Olmanson et al., 2020), suspended sediments (Dekker et al., 2001; Ritchie & Cooper, 1988), water clarity (Olmanson et al., 2008; Topp et al., 2021), and primary productivity (Kuhn et al., 2020). To infer water quality, these studies build models based on relationships between optically active constituent concentrations and their impact on water surface reflectance. These efforts are becoming increasingly accessible due to emerging datasets that match satellite observations with field measurements of water quality parameters for model training and development (Dethier et al., 2020; Ross et al., 2019; Spyrakos et al., 2020), as well as online processing and data storage platforms such as Google Earth Engine (Gorelick et al., 2017).

Here, we present a 36 years analysis of U.S. lake color phenology using LimnoSat‐US, a new analysis‐ready remote sensing data set for inland waters. LimnoSat‐US contains all cloud‐free Landsat observations of U.S. lakes larger than 0.1 km2 between 1984 and 2020. As either a stand‐alone resource, or when combined with existing datasets such as AquaSat (Ross et al., 2019) and RiverSR (Gardner et al., 2021), LimnoSat‐US provides opportunities for novel analyses of remotely sensed, macroscale patterns in U.S. freshwater resources. Through this initial application of LimnoSat‐US, we identify the dominant phenology patterns in U.S. lakes, how those patterns have changed over time, and what lake and landscape level characteristics control the stability of a given lake's seasonal cycle.

Materials and Methods

Database Development

We constructed the LimnoSat‐US database by extracting the surface reflectance values of high confidence water pixels (Jones, 2019) from Landsat 5, 7, and 8 imagery across 56,792 lakes contained within the HydoLAKES database (Messager et al., 2016). While these surface reflectance images were originally developed for terrestrial applications, a growing body of research shows that they can be used to accurately estimate inland water quality parameters and perform on par with water‐specific atmospheric correction algorithms (Griffin et al., 2018; Kuhn et al., 2019; Olmanson et al., 2020). To avoid signal noise from surrounding land pixels and bottom reflectance in shallow waters, we take the median reflectance values from within 120 meters of the deepest point for each waterbody (Shen et al., 2015), where the deepest point refers to the portion of the lake that is furthest away from the lake shoreline. Clouds, cloud shadow, and snow/ice in all images were identified using the Landsat pixel quality flags (Zhu et al., 2015). Observations were removed from the database if any clouds, cloud shadows, snow, or ice were detected within 120 meters of the deepest point, or if the observation contained fewer than nine high confidence water pixels (Jones, 2019).

For each observation within LimnoSat‐US, we calculate water color. Color is an intuitive measure of lake water properties that can be calculated without any knowledge of the inherent optical properties of the water column (Gardner et al., 2021; Giardino et al., 2019; Woerd & Wernand, 2015). For each observation within the LimnoSat‐US database, we quantified water color as both the dominant wavelength (d) and Forel‐Ule index (Wang et al., 2015). Dominant wavelength is calculated by converting reflectance values across the visible spectrum (blue, green, and red) into the color wavelength (nm) as it is perceived by the human eye. By incorporating spectral information across the visible spectrum, this conversion captures nuanced changes in color that may not be evident using a simple metric like the reflectance maximum. The Forel‐Ule Index further breaks this dominant wavelength into 21 discrete color categories that are representative of overall lake typology. A detailed technical description of database development, including atmospheric corrections, sensor standardization, deepest point calculations, and the calculation of dominant wavelength can be found in (Sections S1 and S2).

Seasonal Lake Color Phenology

The development of the LimnoSat‐US database provides novel opportunities for examining macrosystem patterns in U.S. lake dynamics. Clustering analysis is one common approach for extracting patterns from time series datasets that have no a priori assumptions about group membership (Byrnes et al., 2020; Savoy et al., 2019; Warren Liao, 2005). The overall goal of clustering analysis is to partition group membership based on within‐group similarity and between‐group dissimilarity. Here, we apply clustering analysis to time series of lake color to better understand the drivers of variation in lake seasonality over the past 36 years.

Lake color observations generated from the LimnoSat‐US database were filtered to those between May and October to remove missing data caused by snow and ice. Observations were broken into six distinct periods—(1984, 1990], (1990, 1996], (1996, 2002], (2002, 2008], (2008, 2014], (2014, 2020]—and were filtered to those with at least three observations per month per period, resulting in 26,607 lakes with enough data to calculate periodic seasonality for the analyses. The choice of six distinct periods balanced data abundance with temporal specificity (Figure S1). Using six periods results in a median of 10 cloud free observations of lake color per month, per period, for each lake. Within each period, lake color phenology was calculated for both raw dominant wavelength and lake/period z‐normalized dominant wavelength using a Nadaraya‐Watson kernel regression (Nadaraya, 1964; Watson, 1964) implemented with the ksmooth function from the stats package in R (R Core Team, 2019). Application of the kernel regression allowed for the calculation of a weekly color value based on a Gaussian weighted average of all observations within a window of 21 days from the point calculated. Extreme outliers (>4 standard deviations from the lake/period mean) were removed prior to the kernel regression for each series. The resulting time series consist of weekly estimates of lake color from May to October for each lake for each period (Figure 1).

Figure 1

Examples of the calculated seasonal phenologies for three well studied lakes of different trophic states. Landsat observed color points are plotted in the actual color of the observation according to the Forel‐Ule index. Phenologies are composed of one observation per 7 days calculated by taking a Gaussian weighted average of all points ±21 days from each calculated point. Lakes/periods marked by an asterisk were classified as aseasonal and placed in the a priori aseasonal cluster.

Normalization of the time series is critical for accurately clustering lake phenologies using the dynamic time warping (DTW) method described below (Keogh & Kasetty, 2003; Mueen & Keogh, 2016). However, by standardizing the variance across time series, we artificially impose equal seasonal variation between lakes/periods that are relatively monotonic (i.e., aseasonal) and those that show true seasonality in the phenology of their color. Examination of the mean and standard deviation of dominant wavelength for the nonnormalized time series shows that this is particularly problematic for end member lakes on either end of the color spectrum that show very little seasonal variation (σ < 5 nm, Figure S2). This can be seen in Figure 1, where oligotrophic Crater Lake shows minimal seasonality when compared to known eutrophic waterbodies (Lake Mendota and Lake Okeechobee). To address this issue while still following best practices of normalization for clustering analysis, those lakes/periods with a dominant wavelength standard deviation of less than 5 nm were classified a priori as aseasonal. To identify this threshold we examined cut‐offs at the 20th, 30th, and 40th percentiles of standard deviation across lake/period timeseries. The 20th percentile failed to fully capture end‐member lakes, whereas the 40th percentile removed an unnecessary number of time series from the clustering analysis. We therefore chose the 30th percentile (σ = 5.01) as the cutoff threshold to maximize the number of timeseries available for clustering analysis, guarantee that seasonal variation within those timeseries is at least ∼10 nm around the mean color, and effectively classify aseasonal, monotonic, and end‐member lakes into their own grouping.

This process resulted in 109,643 individual time series available for cluster analysis and an additional 46,759 classified a priori as aseasonal. These time series were clustered using DTW (Sakoa & Chiba, 1978) within a partitional clustering framework with barycenter averaging (Sarda‐Espinosa, 2019). DTW allows points within two time series to be compared within a user‐defined window as opposed to using a one‐to‐one comparison found in traditional metrics like Euclidean distance. This elasticity reduces the impacts of noise, minor temporal shifts, and outliers, making it ideal for ecological systems with natural interannual variations (Savoy et al., 2019; Zhang & Hepner, 2017). The final number of clusters was determined by comparing the Davies‐Bouldin (Davies & Bouldin, 1979) and Modified Davies‐Bouldin (Kim & Ramakrishna, 2005) cluster validity indexes (CVI) across iterations ranging from 2 to 8 clusters. The Davies‐Bouldin and Modified Davies Bouldin were chosen because of their computational efficiency and strong performance when compared to other common CVIs (Arbelaitz et al., 2013).

One important validation of clustering analysis is how sensitive final clusters are to sample variations in their input, the idea being that stable, or “universal,” clusters will emerge across differing sampling schemes (Jain & Moreau, 1987). Here, we addressed issues of cluster stability using the Jaccard Similarity Index across 100 iterations of bootstrap sampling of our input time series. At each iteration, the original input time series were sampled with replacement, clustered, and the resulting clustering algorithm used to predict groupings for the original data. The Jaccard Similarity Index was then calculated based on how similar each new cluster was to the corresponding original cluster. The index ranges from 0 to 1, indicating that clusters share all or no members, with values greater than 0.5 generally indicating cluster stability and representativeness of true patterns within the data (Savoy et al., 2019). Significant differences in the distribution characteristics of the final clusters were identified using the non‐parametric Kruskall Wallace Analysis of Variance on Ranks (Hollander & Wolfe, 1973) followed by Dunn's Test with a Bonferroni p‐value correction (Dunn, 1961). These tests provide a conservative (i.e., strict) measure of significant differences and are appropriate for large datasets that don’t meet parametric assumptions of normality.

Finally, we examined the spatial autocorrelation of clusters and the overall stability of individual lake phenologies. Spatial autocorrelation was measured by randomly sampling 30% of the lakes, assigning them their most common cluster, and calculating the proportion of same cluster lakes versus different cluster lakes within 50 km windows moving outward from each lake in the subsample. Lake phenology stability was calculated by examining the number of times a given lake shifted between clusters throughout the six periods of study. Lakes were categorized on a scale from 0 (stable) to 5 (unstable) based on the total number of cluster transitions they made between 1984 and 2020. Lake and landscape level factors from HydroLAKES (Messager et al., 2016) and the Global Lake Area, Climate, and Population database (Meyer et al., 2020) were then used to assess lake characteristics that influence the stability of a lake's seasonal phenology over time. Variables that potentially influence stability were identified through linear regression of lake stability (0–5) on the median value of the lake/climate attribute within each stability class. Those attributes with a coefficient p‐value of less than 0.05 were further examined as correlates with lake stability.

Results

The final LimnoSat‐US database includes reflectance values spanning 36 years for 56,792 lakes across >328,000 Landsat scenes. After initial quality control measures, the database contains over 22 million individual lake observations with an average of 393 ±233 (mean ± standard deviation) observations per lake over the entire study period. While observations date back to 1984, the total number for any given year approximately doubles with the launch of Landsat 7 in 1999 (Figure 2).

Figure 2

Temporal and spatial distributions of satellite observations contained within the LimnoSat‐US database.

Classes of Lake Color Phenology

Our final clustering partitions resulted in one of three membership classes for each lake/period that was not a priori classified as aseasonal (Figure 3). We describe these groups as Spring Greening, Summer Greening, or Bimodal. High mean Jaccard Similarity Indices across bootstrap sampling iterations (0.77, 0.80, and 0.94 respectively) show these clusters are relatively universal, and that regardless of the initial sample, the same lakes are consistently clustered together. Within these clusters, we refer to red‐shifted portions of the time series (increasing values) as greening or green‐shifted and blue shifted portions of the time series (decreasing values) as blueing or blue‐shifted. We highlight this terminology because even though red is the end‐member of the upper wavelengths, the vast majority of the colors do not extend beyond the green portion of the spectrum. It is important to note that while increases in chlorophyll‐a from phytoplankton growth can lead to “greening” patterns, increases in CDOM and non‐algal particles can also cause a lake to move from bluer wavelengths toward the green and red portions of the spectrum. Descriptions of the summary attributes for each cluster are as follows:

Figure 3

Results of cluster analysis for over 26,000 lakes and 156,000 seasonal time series. Black lines represent medians with gray ribbons representing the first–third quartile of each cluster. Clusters are shown both in their (a) z‐normalized form used in the cluster analysis and (b) their raw dominant wavelength form. Distributions of color observations in each cluster (c) are colored by their associated Forel‐Ule Index value and represent the distribution of the human perceived watercolor within each cluster. Note that the range of wavelengths associated with each Forel‐Ule Index value varies.

Spring Greening (n = 55,378, 35.4%): Lake color is green‐shifted in May/June and gradually moves toward the blue end of the spectrum throughout the summer and fall months. Median dominant wavelengths for these phenologies are significantly bluer (p < 0.0001) than those in the Summer Greening, Bimodal, or Aseasonal clusters (median  = 513). They have the highest average coefficient of variation within each individual time series (p < 0.0001), with an average range of 37 nm for a given lake/period compared to 34, 33, and 12 nm for Summer Greening, Bimodal, and Aseasonal clusters, respectively. The distribution of colors within the cluster is concentrated around a mode 498 nm and skewed toward the greener portion of the spectrum.

Summer greening (n = 24,580, 15.7%): Lake color is characterized by gradual greening from May to August after which time it drops toward the blue end of the color spectrum. The distribution of colors shows a mode of 542 nm and a median of 524 nm with a blue‐skewed distribution. On average, each individual time series within this class shows significantly less variation than Spring Greening lakes/periods (p < 0.0001) but no significant difference from Bimodal lakes/periods.

Bimodal (n = 29,685, 19.0%): Lake color is most green‐shifted in May/June and again in September/October, with a somewhat blue‐shifted phase in the intervening months. Phenologies within this cluster are significantly more green‐shifted (p < 0.0001) than lakes within either the Spring or Summer Greening clusters and show less variation (p < 0.0001) than those in the Spring Greening clusters. The distribution of colors is concentrated around 553 nm with a much less pronounced peak at 507 nm.

Aseasonal (n = 46,759, 29.9%): The overall color distribution of this cluster is distinctly bimodal, with a primary mode at 559 nm and a secondary mode at 492 nm. This bimodal distribution, combined with the small variance in any given lake/period in the cluster, suggests it contains predominantly blue and predominantly green time series with very few observations in the intermediate green/blue space common within the three other clusters. The cluster also contains both the most green‐shifted and most blue‐shifted time series included within the analyses. Because of the crisp partition contained within the cluster and the ecological significance of blue versus green aseasonal time series, we further partition this cluster into Aseasonal (Blue) (n = 15,934) and Aseasonal (Green) (n = 30,825) lakes for the remainder of the analysis. Time series with a median dominant wavelength less than or greater than the anti‐mode of the distribution (525 nm) are considered Aseasonal (Blue) and Aseasonal (Green) respectively.

Lake Stability Over Time

Aseasonal Green lakes showed the most stability over time, with an average of 73% ± 6% (mean ± standard deviation) of lakes remaining within the cluster between consecutive time periods. Aseasonal (Blue) and Spring Greening clusters showed similar retention rates of 57% ± 17% and 57% ± 9% respectively, while Bimodal and Summer Greening showed similar retention rates of 46% ± 8% and 45% ± 7% (Figure 4). However, of these, only the differences between Aseasonal (Green) and Bimodal/Summer Greening clusters were statistically significant at a 95% confidence interval. For Spring Greening, Aseasonal (Green), and Aseasonal (Blue) distributions, the number of lakes retained between each period was significantly higher than the number of lakes that transitioned to a different cluster (p = 0.047, p = 0.007, and p = 0.0001 respectively). Summer Greening and Bimodal clusters showed no significant difference between the proportion of lakes retained and lakes that transitioned to other clusters, indicating less stability than the other three classes. However, these transitions showed distinct patterns, with lakes transitioning more commonly between similar clusters. As an example, on average 27% of Summer Greening lakes transitioned to Spring Greening lakes between periods, but only 4% of Summer Greening lakes transitioned to Aseasonal (Green). Similarly, less than 0.2% of lakes in Aseasonal (Green) and Aseasonal (Blue) transitioned between the two clusters in any two consecutive periods, indicating that state shifts between dominantly blue lakes and dominantly green lakes are very uncommon.

Figure 4

Sankey diagram showing the distribution of lake phenology transitions between periods. Each ribbon is proportional to the number of lakes that moved from one cluster class to another.

Lake stability, or the number of times a lake moved from one cluster to another (ranging from 0 transitions to 5), showed that lakes with three transitions were most common (n = 6,458) and lakes with five transitions least common (n = 1,254) (Figure S3). We also calculated the number of unique clusters a lake occupied throughout its transitions. For instance, a lake could change states between all five periods, giving it a stability score of five, but only be changing between two of the potential five clusters, giving it two unique states. Of the 26,067 lakes, 4,339 (16.6%) remained within the same cluster through all periods while only 21 (<0.1%) occupied all five clusters at some point. For those lakes in between, lakes occupying two distinct states (n = 11,091; 42.5%) were most common followed by three states (n = 8,942; 34.3%) and four states (n = 1,674; 6.5%) respectively. Linear regressions between lake and landscape level metrics with overall lake stability showed significant relationships (p < 0.01) with 5 out of 26 possible metrics (Table S1), although some of these metrics have significant cross‐correlation (Figure 5).

Figure 5

Correlation matrix for variables from HydroLakes and GLPC used to analyze drivers of lake phenology stability. Both the size and color of the squares indicate the strength of the Pearson's correlation coefficient between any two given variables.

Discussion

Lake Seasonal Phenology Types

Existing paradigms regarding the seasonality of lake color are generally derived from individual lakes with rich sampling histories of water quality observations; however, these long‐term field records are rare and limited to a small subsample of lakes (Stanley et al., 2019). While these data‐rich study lakes are essential for understanding fine‐scale ecosystem processes, they lack the spatial coverage to generalize across entire landscapes (Collins et al., 2019; Soranno et al., 2014). Within our clustering analysis, we found that lake color phenology can largely be categorized as Aseasonal, Spring Greening, Summer Greening, or Bimodal. These phenologies show distinct regional patterns and spatial auto‐correlation, with the probability of two lakes being in the same cluster showing a significant relationship to the distance between those two lakes (p < 0.0001) up to a distance of ∼1,500 km (Figure 6b).

Figure 6

(a) The modal cluster within each 100 × 100 km grid across time periods. Mixed grids are those where there is no dominant cluster (i.e., two or more clusters are equally prevalent). (b) The frequency of same cluster pairs to different cluster pairs using each lake's modal cluster. The frequency distributions were calculated within 50 km windows for a random sample of 30% of the study lakes. The dotted line represents the expected frequency if the distribution was random without any spatial autocorrelation.

Each cluster has a unique distribution of dominant wavelengths (Figure 3), which suggests that the timing of seasonal variation in color is connected with lake biogeochemistry. This conclusion is supported by long‐standing models of freshwater phytoplankton succession (Sommer et al., 1986) and observations of annual cycles of chlorophyll‐a, a proxy for phytoplankton biomass (Winder & Cloern, 2010). Oligotrophic temperate lakes often show the archetypal pattern of a spring phytoplankton bloom followed by low summer concentrations. This was the dominant phenology in our observations (35.4%), which is in‐line with a study of 125 aquatic systems that found that nearly half of the sites displayed a dominant 12‐months cycle with one phytoplankton peak per year (Winder & Cloern, 2010). As nutrient availability increases, eutrophic lakes tend to experience discrete phytoplankton blooms in the spring and late‐summer/fall (Marshall & Peters, 1989). This pattern is captured in our Bimodal cluster, where the raw dominant wavelength values are significantly greener than those in any other cluster except for Aseasonal (Green). The summer‐greening cluster captures eutrophic to hyper‐eutrophic lakes featuring prolonged summer blooms with highly variable summer algal concentrations (Carpenter et al., 2020; Huisman et al., 2018). The characterization of Bimodal and Summer Greening lakes/periods as eutrophic is further supported by the low levels of variation we observe in dominant wavelengths when compared to Spring Greening lakes/periods. Dominant wavelength saturates with high amounts of suspended matter, chl‐a, and/or CDOM (Bukata et al., 1997), meaning that highly productive, algae‐filled lakes with significant amounts of these constituents would show low variation as dominant wavelength saturates. It is also possible that some lakes in these categories are dystrophic CDOM‐dominated lakes, as they include some of the more red‐shifted (brown) waterbodies within the study.

The proportion of lakes that fall within different clusters does not show an overall trend over time; however, since the 1996–2002 period, the number of lakes classified as either Bimodal or Aseasonal (Blue) have increased while the number classified as Spring Greening have been decreasing (Figures 4 and 6). Much of the increase in Aseasonal (Blue) lakes is concentrated in the Pacific Northwest and occurred prior to 2008, whereas the decrease in Spring Greening Lakes has predominantly occurred in higher‐latitude lakes that may be more sensitive to changes in snowmelt and ice cover regimes which control nutrient and sediment fluxes that influence lake productivity (Gerten & Adrian, 2002; Sharma et al., 2019). Patterns in the Aseasonal (Green) cluster show much less variation both spatially and temporally, being largely concentrated in the agriculturally dominated central and northern plains and showing no distinct temporal pattern in quantity. While the increase in Aseasonal (Blue) lakes is potentially indicative of reduced sediment and nutrient inputs in certain parts of the country, the increase in Bimodal lakes, when taken with its close match to eutrophic phytoplankton succussion patterns, indicates increases in lake productivity across portions of the U.S. since the mid 1990s. This pattern supports recent research showing a transition from bluer lakes to murky chlorophyll‐a and CDOM‐dominated lakes throughout the US between 2007 and 2012 (Leech et al., 2018), as well as work highlighting that while overall clarity in U.S. lakes has increased since 1984, these increases largely took place prior to the mid‐1990s, after which time clarity improvements began to plateau (Topp et al., 2021). However, dominant wavelength, and optical water color more generally, is controlled by a variety of optically active water color constituents in addition to phytoplankton (Mobley, 1994), and partitioning these optical components is beyond the scope of this analysis. The result does, however, merit further research using a database like LimnoSat‐US to examine countrywide trends in lake chlorophyll‐a content.

Factors Influencing Lake Stability Over Time

Lake stability, or the number of times a lake moved between clusters during the study period, showed significant relationships with multiple lake and landscape level metrics from HydroLAKES and the GLCP database (Figure 7, Table S1). These relationships can generally be categorized as either hydrological properties or landscape properties. Important hydrological properties related to stability include lake size and lake inflow (both positively correlated with stability). This result supports existing research suggesting that larger water bodies are less reactive to perturbations than smaller, shallower lakes that can fluctuate among multiple productivity regimes (Scheffer & van Nes, 2007). We also find that hydrologically dynamic lakes are consistently less stable, with lakes showing large interannual variations in seasonal surface extent exhibiting less stability. It is likely that these hydrologically dynamic lakes are more sensitive to seasonal variations in runoff and resuspension of lakebed sediments leading to large interannual variations in nutrient and sediment load.

Figure 7

Lake and landscape level metrics that showed the most significant relationships with stability, or the number of times a given lake moved from one cluster to another between periods (p < 0.01 with the exception of lake inflow, p = 0.019). Center bars represent median values while boxes span the first‐third quartiles.

The landscape‐level metrics that showed the strongest relationship with lake stability were catchment population and elevation (p < 0.01) followed by mean temperature and mean monthly precipitation (p < 0.05). Similarly, for the subset of these variables where we had observations at annual timescales, we found that coefficients of variation between years (interannual variation) of these metrics showed strong linear relationships to stability. The impact of these landscape‐level metrics on stability supports work showing that lakes integrate surrounding climatic and land cover changes (Rose et al., 2017). These results are of particular interest for relatively pristine high‐elevation lakes that will be disproportionately impacted by changing precipitation and temperature regimes through climate change (Oleksy, Baron, et al., 2020; Oleksy, Beck, et al., 2020). Finally, we found that lakes in catchments with higher populations were generally more stable; however, lakes in catchments with high variation in population (likely increasing urban areas) showed less stability. Overall, our examination of landscape level metrics shows that the stability of a lake often follows the stability of its environment, with lakes subject to interannual variations in climate or anthropogenic stressors generally showing less stability in their overall seasonal phenology.

Conclusion

Remote sensing has the capability to substantially increase macroscale understanding of aquatic ecosystem processes. Here, we contribute to a growing body of inland water remote sensing resources with LimnoSat‐US, which contains >22,000,000 remotely sensed lake observations. Prior to this study, large‐scale analyses of lake phenologies were limited to dozens to hundreds of waterbodies (Ho et al., 2019; Marshall & Peters, 1989; Winder & Cloern, 2010). Here, we were able to analyze U.S. summer lake color phenology across more than 26,000 lakes over 36 years, showing both temporal and spatial patterns and trends, as well as linking phenology to lake and landscape‐level metrics. Better understanding the full distribution of lake phenology will allow for more accurate scaling of global nutrient and carbon cycling. While the analysis presented here relies simply on lake color, combining LimnoSat‐US with databases such as AquaSat (Ross et al., 2019), RiverSR (Gardner et al., 2021), and LIMNADES (Spyrakos et al., 2020), will allow for more explicit modeling and analysis of specific water quality components, allowing researchers to partition the patterns observed here into optically active water quality components including chlorophyll‐a, suspended sediments, and CDOM.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Supporting Information S1

Click here for additional data file.