Long-term trends and drivers of water color in Missouri reservoirs

Contrasting water quality trends are occurring within and across North America, with waterbodies experiencing increasing phytoplankton blooms, increasing dissolved organic matter, or both, while others are becoming clearer and bluer; dramatically changing water color. To assess the spatial and temporal variability in water color, we quantified trends in satellite-derived dominant wavelength (λ_d) from 1984 to 2020 for 484 reservoirs across the state of Missouri using the LimnoSat-US dataset. Currently, the vast majority of Missouri reservoirs are classified as green and within a range (538–555 nm) that lies closer to the brown color endmember. Nearly one-third of reservoirs (n = 159) experienced significant temporal shifts in water color, with more (n = 91) negative (e.g., bluer) than positive (n = 68) λ_d trends. Linear mixed-effect models indicate that periods of extreme wetness and drought are associated with browner and bluer waters, respectively, and boosted regression trees further reveal that waterbody and watershed characteristics are important predictors for water color trends. We also analyzed trends in summer water quality (WQ) parameters from two long-term monitoring programs to evaluate independent and synchronous changes with λ_d. We provide analyses showing that particulate inorganic matter and Secchi depth most strongly correlated with λ_d, and total nitrogen and total phosphorus concentrations that are not typically associated with satellite-derived data have greater co-variance with λ_d than chlorophyll a. Although bluer waters often reflect reductions in inorganic turbidity, our findings show that they can at times coincide with increases in chlorophyll a. We further demonstrate that while λ_d trends broadly align with changes in water quality, co-occurring water quality and color trends in Missouri reservoirs at times defy a simplistic canonical interpretation, particularly in eutrophic waterbodies where changes in nutrient concentrations, chlorophyll a and water color can occur independent of each other. Our results help explain some of the previously observed heterogeneous controls on water color and emphasize the importance of integrating water quality data alongside commonly used landscape and morphological features.

1 Introduction

Contrasting water quality trends are occurring within and across North America, with waterbodies experiencing oligotrophication (Sillen et al., 2024), eutrophication (Schindler et al., 2012), brownification (Monteith et al., 2007; Roulet and Moore, 2006), or both simultaneously (Leech et al., 2018). Water color, as perceived by the human eye, is one of the oldest indicators of water quality parameters and is closely linked to productivity and trophic state of aquatic ecosystems (Topp et al., 2021). A number of commonly measured water quality variables can be estimated directly from direct or remotely-sensed measures of color including chlorophyll a (chl a; Smith et al., 2020), total suspended solids (TSS; Ondrusek et al., 2012), and colored dissolved organic matter (CDOM; Cao et al., 2018). Indeed, much of our understanding of canonical shifts in water quality are intertwined with changes in color. For instance, high nutrient loads into waterbodies from urban, agricultural, and industrial sources result in increased chl a (proxy for phytoplankton biomass) and greener waters (Dodds et al., 2009). Increases in CDOM due to changes in climate, hydrology, land cover, and atmospheric deposition (Erlandsson et al., 2008; Freeman et al., 2001; Monteith et al., 2007) have led to brown waters (Leech et al., 2018). Conversely, factors such as lake acidification (Charifson et al., 2015), low precipitation (Hongve et al., 2004), increased rates of filter-feeding by zebra mussels (Binding et al., 2007), or reductions in nutrient loading (Jeppesen et al., 2005) lead to increased water transparency and bluer waters. The nutrient-color paradigm (Williamson et al., 1999) conceptualizes the interaction between total phosphorus (TP) and CDOM in lakes where blue waters are oligotrophic (low TP and CDOM), green waters are eutrophic (high TP, low CDOM), and brown waters are either dystrophic (low TP, high CDOM) or mixotrophic (‘murky’, high TP and CDOM; Oleksy et al., 2024; Williamson et al., 1999).

Water color is strongly representative and accessible through remote sensing making it an attractive indicator for monitoring changes in the ecological state and environmental conditions of aquatic ecosystems (Shen et al., 2025). Satellite-based monitoring offers a cost-effective alternative, providing broad spatial and temporal coverage, but its application to inland waters has been hindered by the limited capabilities of satellite sensors not specifically designed for coastal and inland waters applications and by challenges such as atmospheric correction, land adjacency, and bottom reflectance (Mouw et al., 2015). Recent advances, including improved atmospheric correction algorithms (Castagna and Vanhellemont, 2025), machine learning approaches (Smith et al., 2020), and robust analysis-ready databases (King et al., 2024; Topp et al., 2020), collectively aim to improve data quality and accessibility. Satellite-derived dominant wavelength (λ_d) has emerged as a widely used descriptor because it captures the perceived hue of water, is readily derived from virtually all earth observing satellites, and provides a consistent metric for studying long-term water quality changes at broad scales (Gardner et al., 2021; Yang et al., 2022). Studies using λ_d have revealed the importance of climate, landscape characteristics and waterbody morphology in driving water color trends (Cao et al., 2023; Liu et al., 2024; Oleksy et al., 2022; Shen et al., 2025).

Despite large scale efforts (Gardner et al., 2021; Yang et al., 2022), the factors driving changes in water quality and color may differ both among individual waterbodies and across geographic regions (Cao et al., 2023). As a result, more focused regional scale studies may be more effective for describing parallel changes in water quality and water color. However, the scarcity of long-term in situ data has limited our ability to assess these dynamics over broad spatial and temporal scales. This study integrated 37 years of satellite-derived water color observations from the analysis-ready LimnoSat dataset with detailed field-based measurements from long-term water quality monitoring programs in the state of Missouri to 1. quantify secular and synchronous changes in summer water color and water quality in Missouri reservoirs over the past several decades and 2. identify climate, landscape characteristics, and waterbody morphology factors driving the spatial and temporal trends. Inspired by the nutrient-color paradigm, this study examined covariation not only between optically active parameters (e.g., Chl a, TSS) and water color, but also tests the suitability of λ_d as a predictor of long-term changes in total phosphorus (TP) and total nitrogen (TN). Previous studies have demonstrated that non-optically active components such as nutrients can be retrieved from satellite-derived surface reflectance (Gao et al., 2015; Guo et al., 2024; Li et al., 2017; Sun et al., 2014; Xiong et al., 2022), and that water color can be used as a metric in trophic state classification (Meyer et al., 2024; Webster et al., 2008). However, to our knowledge, this is the first study that examines how long-term, satellite-derived water color responds to multi-decadal changes in nutrient concentrations across hundreds of reservoirs. By compiling a uniform dataset from open data repositories and statewide long-term water quality programs, our study not only aims to unravel the key drivers of regional trends in water color and quality in the study region, but also demonstrates the potential for integrating satellite-based measurements with traditional monitoring to strengthen future research; particularly in times of reduced funding for environmental monitoring programs.

2 Materials and methods

2.1 Study area and geospatial analysis

Most of Missouri’s (MO) waterbodies are human-constructed reservoirs embedded across the state’s six Environmental Protection Agency (EPA) ecoregions (Ozark Highlands, Central Irregular Plains, Western Corn Belt Plains, Interior River Valleys and Hills, the Mississippi Alluvial Plain, and Mississippi Valley Loess Plains; Figure 1A). Geospatial data for individual reservoirs, published by the MO Department of Natural Resources (MDNR) was separately matched in space with the in situ water quality monitoring datasets and the remote sensing dataset described below. The MDNR identifies a total of 2,547 reservoirs, of which 484 have remotely sensed data with ≥ 20 consecutive years of summer color observations (Section 2.2); 298 have water quality data (Section 2.3), and 160 have co-located water quality and remotely sensed time series (Figure 1A).

Figure 1

Figure 1. Map of (A) Missouri (United States) reservoirs showing available in situ and/or LimnoSat time series embedded across the state’s six Environmental Protection Agency (EPA) ecoregions (Ozark Highlands, Central Irregular Plains, Western Corn Belt Plains, Interior River Valleys and Hills, Mississippi Alluvial Plain, and the Mississippi Valley Loess Plains), and (B) predominant watershed land cover types of the studied reservoirs based on 2021 NLCD data.

Watersheds were delineated using the Python package PySheds and USGS 60 m digital elevation data. Watershed land cover was extracted from the 2021 release of the National Land Cover Database (NLCD), developed primarily by the U.S. Geological Survey (USGS) in collaboration with the Multi-Resolution Land Characteristics (MRLC) Consortium. Watershed land cover was treated as a static variable in this study as a previous state-wide study demonstrated <2% change in any land use category over 15-year period that overlaps with this study (Bhattacharya et al., 2022). Geomorphological data, specifically elevation, area, perimeter, and elongation ratios were derived for each reservoir and watershed. The elongation ratio is defined as the ratio of the diameter of a circle with an equivalent area to that of the reservoir/watershed to its maximum length (i.e., the maximum distance between all coordinate pairs along the perimeter; Sukristiyanti et al., 2018). High elongation ratios therefore represent more circular reservoirs/watersheds and generally coincide with flat land and low relief. Smaller elongation ratios correspond to more narrow and dendritic reservoirs/watersheds with high relief and steep slopes (Sukristiyanti et al., 2018). All geospatial analysis was performed in Python, and the functions and libraries used for each calculation are listed (Supplementary Table S1).

Hydrometeorological variations in space and time were analyzed using the Palmer drought severity index (PDSI). The PDSI index quantifies the severity of wet and dry conditions using estimates of relative soil moisture conditions based on temperature, precipitation, and evapotranspiration anomalies. PDSI values less than −3 indicate severe drought conditions, and PDSI values greater than 3 indicate very moist conditions. PDSI for summers from 1990 to 2020 were obtained from the daily high-spatial resolution Gridded Surface Meteorological (gridMET) dataset available on the Google Earth Engine (GEE) platform.

2.2 Water color

Dominant wavelength (λ_d), a metric derived from standardized surface reflectance data, was obtained for 484 waterbodies in Missouri from the LimnoSat-US database (Topp et al., 2020). The main criterion for selecting the studied reservoirs was the availability of 20 consecutive years of data between 1984 and 2020. For inland waterbodies in the U.S. larger than 0.1 km², LimnoSat contains all cloud-free Landsat (5, 7, and 8) surface reflectance (T1-SR) data where for each waterbody and timestamp, the median reflectance is derived from pixels within 120 m of the Chebyshev Center (i.e., the deepest point that is furthest away from the shoreline). Landsat 5 and 7 imagery have been atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS), while Landsat 8 imagery has been corrected using the Landsat Surface Reflectance Code (LaSRC). Bands for each sensor were standardized across time and between satellites (Topp et al., 2021), ensuring the reliability of long-term trend analyses (Shen et al., 2025). Dominant wavelength is calculated by converting surface reflectance values into the chromaticity color space, following Wang et al. (2014). For each reservoir and year, the mean λ_d was computed across the stratified season (May 1 to October 31) to align with the water quality sampling period. The number of valid summer λ_d observations per reservoir-year is provided in the Supplementary Material (Supplementary Figure S1).

A variety of λ_d thresholds have been used in the literature to delineate blue from green waterbodies and green from brown waterbodies (Oleksy et al., 2022). In the Missouri LimnoSat data, blue reflectance is maximal where λ_d ≤ 495 nm (blue reservoirs), red reflectance is maximal where λ_d ≥ 575 nm (brown reservoirs), such that green reservoirs have λ_d between 495 nm < λ_d < 575 nm (Supplementary Figure S2). Within the green spectral range, blue reflectance steadily declines while red reflectance steadily increases (Supplementary Figure S2), such that decreasing trends in λ_d can be described as blueing and increasing trends in λ_d can be described as browning.

2.3 In situ water quality dataset

Water quality parameters for 298 waterbodies in Missouri from 1984 to 2020 were obtained through the Statewide Lake Assessment Program (SLAP) and the community science led Lakes of Missouri Volunteer Program (LMVP), aligning with the LimnoSat-US data availability. In each waterbody, surface water quality samples are collected alongside vertical profiles on 3 or 4 occasions between May and October each year, in close proximity to the deepest part of the reservoir (i.e., near-dam locations). Water temperature and dissolved oxygen profiles, recorded in meter increments (1984–2016) or continuous profiles (2017–2020) from the surface to the bottom, were obtained using multiple YSI sondes over the years (see Supplementary Table S2). These data were used to calculate average water temperature and dissolved oxygen concentrations in the epilimnion (above the top of the metalimnion) and the hypolimnion (below the bottom of the metalimnion). The thermal stability index, Brunt-Väisälä buoyancy frequency (N²; s^-2), was calculated based on the density gradient of the water column in the R statistical environment, version 4.4.2 R Core Team, (2024); (Supplementary Table S1). Thermocline depth at 0.3 kg m^-3/m density threshold was derived from water temperature profiles. Water transparency (Secchi) was assessed using a Secchi disk. Surface water samples refer to either samples from 0.5 m depth or epilimnetic integrated samples (Supplementary Figure S3). Water samples were stored in high density polyethylene (HDPE) bottles, placed into coolers, processed, and then frozen until analysis. After freezing, they were analyzed by standard methods (Supplementary Table S2) for particulate organic matter (POM; mg/L), particulate inorganic matter (PIM; mg/L), total suspended solids (TSS; mg/L), total phosphorus (TP; µmol/L), total nitrogen (TN; µmol/L), uncorrected chlorophyll a (Chl a; µg/L), and dissolved organic carbon (DOC; µmol/L). Sampling and analytical methodology were consistent throughout, except for chl a, TP, and TN (Supplementary Table S2). We also determined the ratios of TN and TP (TN/TP), POM and TSS (POM/TSS), PIM and TSS (PIM/TSS). A statistical summary for all in situ variables is provided in the Supplementary Material (Supplementary Tables S3–S5).

2.4 Data analysis

2.4.1 Trend analysis

All water quality and λ_d trend analyses were performed on reservoir specific annual data averaged across the thermally stratified season (May 1 to October 31). Rates of change were estimated using Sen’s slope (Q; see Supplementary Table S6 for statistical references), and the statistical significance of resultant trends was assessed using the non-parametric Mann-Kendall trend test (MK; Supplementary Table S6). Secular trends of λ_d (i.e., Δλ_d) are limited to reservoirs with a minimum of 20 consecutive years of data between 1984 and 2020 (n = 484, Supplementary Figure S4). Missing values within the time series were imputed using linear interpolation, allowing a maximum gap of two consecutive missing years. No interpolation was performed for missing values at the start or end of a time series. To assess differences in the most recent λ_d and Δλ_d across ecoregions and watershed land cover types, we used nonparametric Kruskal-Wallis tests. When significant differences were found, pairwise Dunn’s post hoc tests with Bonferroni-adjusted p-values were applied to identify specific group differences.

Secular trends of water quality parameters were derived as above but restricted to reservoirs with at least 10 consecutive years of data between 1984 and 2020 (Figure 2). Not all reservoirs (n = 298) have the required time series available for every variable (Supplementary Figure S4). To better capture the complexity of water quality dynamics in Missouri reservoirs, two additional analyses were performed. First, linear regression of resultant water quality slopes, regardless of significance, was conducted in pairwise fashion using parameters with >100 co-occurring time series (TP, TN, Chl a, TSS, PIM, POM, Secchi disk). Second, a frequency analysis is presented that counts how often significant water quality trends co-occur within a reservoir, including cases of contrasting trends (e.g., increasing TN and decreasing TP).

Figure 2

Figure 2. Flowchart illustrating the processing and analysis of climatic (Palmer Drought Severity Index; PDSI), remotely-sensed (LimnoSat), in situ water quality (WQ), reservoir morphology and watershed land cover data. LMM refers to linear mixed-effect models, BRT refers to boosted regression tree, and n refers to the number of reservoirs included in each dataset or analysis. Water quality parameters include chlorophyll a (chl a), Secchi disk depth (Secchi), total nitrogen (TN), total phosphorus (TP), total nitrogen to total phosphorus ratio (TN/TP), particulate organic matter (POM), particulate inorganic matter (PIM), total suspended solids (TSS), surface water temperature (wtemp), hypolimnion water temperature (hypo wtemp), buoyancy frequency (buoy freq), thermocline depth (thermocline), dissolved oxygen (DO) profile metrics (epilimnion DO, hypolimnion DO, oxycline and hypoxycline), and dissolved organic carbon (DOC).

2.4.2 Remotely sensed and in situ water quality matchup analysis

To evaluate how well remotely-sensed water color serves as a proxy for water quality and trends, two analyses are performed using a subset of water quality variables that serve as trophic status indicators (Chl a, TN, TP, TSS, POM, PIM, Secchi, DOC). First, λ_d is compared to collocated water quality observation using a 1-day time window from satellite overpass. The final dataset consists of 1704 λ_d and near-concurrent in situ water quality parameters across multiple reservoirs. Pairwise Spearman rank correlations were used to evaluate associations between λ_d values and water quality parameters. To assess differences across wavelengths, water quality data were grouped by λ_d values binned at 10-nm and tested for statistically significant differences using nonparametric Kruskal-Wallis tests. When significant differences are detected, pairwise Dunn’s post hoc tests with Bonferroni-adjusted p-values are applied. Second, parallel trends in λ_d and water quality variables are examined in a subset of reservoirs (n = 22) with significant λ_d trends and sufficiently long (≥10 years) coincidental water quality data. To assess the degree of correspondence between changes in λ_d and water quality across reservoirs and water quality variables, a simple alignment score was developed as follows. For each reservoir and water quality parameter, a score of +2 is assigned when both trends are statistically significant and in the same direction (e.g., decreasing λ_d and decreasing TP). A score of +1 is assigned when the water quality trend is non-significant but in the same direction as Δλ_d, thus the alignment score places double the numerical emphasis on commensurate trends that are statistically significant (p > 0.05) versus those that are not. Scores of −1 and −2 are assigned when water quality are in the opposite direction of Δλ_d (e.g., decreasing λ_d and increasing TP) and are not significant and significant, respectively. Resultant scores were summed across reservoirs and variables and normalized to a scale of −1, 1. Therefore note that the alignment score does not summarize magnitudes or overall direction of trends, but rather whether or not resultant follow our canonical interpretation of shifting water quality metrics on color. Note that Secchi depth trends are interpreted inversely, with increases indicating improved water transparency.

2.4.3 Boosted regression trees

We employed Boosted Regression Trees (BRTs) to quantify the relative contributions of static spatial controls, geomorphological (waterbody elevation, waterbody and watershed area, waterbody and watershed elongation ratio, Warea/Rarea) and watershed land cover (percent cover of urban, forest, and agriculture lands in the watershed) to variation in dominant wavelength slope (Δλ_d; Sen’s slope rate of change). Our final dataset included 160 reservoirs with co-located LimnoSat, geomorphological and land cover data, sampled between 1984 and 2020 (Figure 2). The BRT method is a machine learning technique that combines regression decision trees and a boosting algorithm (Elith et al., 2008). At each step, a new tree that best reduces the loss function is fitted to the residuals without changing the existing trees as the model is enlarged (i.e., stagewise). This approach eliminates the need for prior data transformation or outlier elimination, fits complex nonlinear relationships between explanatory and response variables, handles interaction effects between predictors, and accommodates both categorical and numeric variables, as well as missing data (Elith et al., 2008). BRTs were implemented in the R statistical environment, version 4.4.2 R Core Team, (2024) using the gbm.step function, which incorporates cross-validation to identify the optimal number of trees and implements the procedures described by Elith et al. (2008). The functions and libraries used are listed in the Supplementary Material (Supplementary Table S1).

In BRT modelling, three hyperparameters need to be defined and tuned: The learning rate (lr), which determines the contribution of each tree to the final fitted model; tree complexity (tc), which controls the size of trees and whether interactions between variables should be considered; and bag fraction (bf), which specifies the proportion of data, without replacement, from the full training set to be selected at each step. The number of trees (nt) is determined from the combination of lr and tc and, in practice, smaller learning rates and larger tree complexities increase the number of trees. To select the appropriate structure of the BRTs, we fitted 30 combinations with varying values for lr (0.05, 0.01, 0.005, 0.001, 0.0005), tc (1, 5, 10) and bf (0.5, 0.75). Model potential for overfitting and performance was assessed using an 80:20 train:test data split with a 10-fold cross-validation method. First, models with the lowest prediction deviation (i.e., deviance between observed values from the training set and predicted values from the test set) were selected. Then, the performance of these models was further evaluated using three metrics: mean absolute error (MAE; Equation 1), root mean square difference (RMSE; Equation 2), and coefficient of determination (R²; Equation 3). Final model selection was based on the combination of lowest MAE and RMSE and highest R² value (i.e., higher R² reflects better predictive performance). The metrics were calculated as follows:

$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n}\left|P_i-O_i\right|(1)$

$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\left(P_i-O_i\right)}^2}(2)$

$R^2=\frac{{\displaystyle {\sum }_{i=1}^n}{\left(P_i-\overline{O}\right)}^2}{{\displaystyle {\sum }_{i=1}^n}{\left(O_i-\overline{O}\right)}^2}(3)$

where n is the number of observations, $P_i$ and $O_i$ are the predicted and observed values at i, respectively, and $\overline{O}$ is the mean of the observed data. MAE indicates how close the prediction is to the observed values, while RMSE represents the standard deviation of the residuals. R² represents the predictive accuracy of models; it approaches 1 when prediction becomes optimal. The final optimal values of lr, tc, and bf were set to 0.0005, 10, and 0.75, respectively, with a maximum number of trees = 3,100 and a Gaussian error distribution (Supplementary Table S7).

The selected model was used to calculate the mean relative importance (RI) and cumulative contribution of each predictor in explaining the variability of Δλ_d. The RI is quantified by analyzing the number of times a predictor is selected for splits in regression trees and the improvement it brings to the model’s predictive accuracy. The RI of each variable is given by dividing its importance by the sum of importance values across all explanatory variables, collectively summing to 100% (Elith et al., 2008). Higher RI values indicate a stronger influence on the response variable. The cumulative contribution reflects the combined effect of multiple predictors in explaining Δλ_d. It is calculated by ranking predictors in descending order of their mean RI and computing their cumulative contributions as a percentage of the total importance. Significant predictors based on permutation tests of variable importance were used to calculate Shapley values, which quantify how each predictor contributes to predictions across all observations. Finally, SHAP-based dependence plots were then used to interpret and visualize how each predictor’s contribution varies across its observed range. Positive SHAP values indicate a contribution to increasing Δλ_d, while negative values indicate a contribution to decreasing Δλ_d across the study reservoirs. The functions and libraries used are listed in the Supplementary Material (Supplementary Table S1).

2.4.4 Linear mixed-effect models

To quantify the influence of hydroclimate (e.g., drought and wetter climate) on water color (λ_d), we employed linear mixed-effect models (LMMs; Figure 2; see Supplementary Table S6 for statistical references). Since hydroclimatic drivers vary through time, LMMs provide an appropriate framework to model their temporal effects while accounting for repeated observations within reservoirs. In all models, λ_d was the response variable, with reservoir identity included as a random intercept to account for inter-reservoir variability in water color. Fixed effects included PDSI, ecoregion, and their interaction. The inclusion of the interaction term allows us to account for climatic effects that might differ by ecoregion. To capture possible underlying temporal trends in λ_d across reservoirs, year was included as a numeric covariate. The variance explained by LMMs was based on marginal and conditional adjusted R² (Nakagawa and Schielzeth, 2013). Marginal adjusted R² (R²m) represents the variance explained by a fixed term, and conditional R² (R²c) represents the variance explained by both fixed and random terms. A visual examination of diagnostic plots was applied to determine the model’s goodness of fit. p-values and 95% confidence intervals were calculated using Wald t-tests with Satterthwaite’s approximation of degrees of freedom. All models were implemented in the R statistical environment, version 4.4.2 R Core Team, (2024). The functions and libraries used are listed in the Supplementary Material (Supplementary Table S1).

3 Results

3.1 Color and trends in Missouri reservoirs

The majority of Missouri reservoirs are green (n = 458, 94%); brown reservoirs (λ_d > 575 nm, n = 24, 5%) are more common than blue reservoirs (λ_d < 495 nm, n = 3, 1%). Reservoir color varies by ecoregion and land use (Figure 3). Reservoirs in the largely forested Ozark Highlands have λ_d values that are bluer than all other ecoregions (Kruskal-Wallis; H = 72.58, p < 0.01; Figure 3B), and reservoirs surrounded by agriculture and wetland watersheds have λ_d values that are browner than reservoirs in forested, mixed, and urban watersheds (Kruskal-Wallis; H = 94.03, p < 0.01; Figure 3D).

Figure 3

Figure 3. (A) Distribution of reservoir dominant wavelength (λ_d) across the state’s six Environmental Protection Agency (EPA) ecoregions, box and whisker plot of (B) λ_d in recent years (i.e., last year for each reservoir) by ecoregions, (C) dominant wavelength slope (i.e., change; Δλ_d), by ecoregion, (D) λ_d in recent years by watershed land cover types, and (E) Δλ_d by watershed land cover types. Different letters represent the significant differences (p-values ≤ 0.05) in λ_d and Δλ_d between ecoregions (B,C) and watershed land cover types (D,E) based on pairwise Dunn’s tests with Bonferroni-adjusted p-values.

Nearly one-third (n = 159) of Missouri reservoirs had statistically significant shifts in reservoir color (Δλ_d, MK analysis, p < 0.05). Amongst these reservoirs, negative Δλ_d are approximately 15% (n = 91) more common than positive Δλ_d (Table 1). Significant shifts in reservoir color, regardless of direction, are almost exclusively confined to shifts within the green spectrum (Figure 4A). Negative Δλ_d reflects shifts from the browner end toward greener hues, whereas positive Δλ_d reflects shifts from the bluer end toward greener hues, rather than transitions toward blue or brown. The number of reservoirs exhibiting negative trends exceeded those exhibiting positive trends in three of Missouri’s six ecoregions (Table 1); however, Δλ_d did not co-vary with ecoregion (Kruskal-Wallis; H = 2.49, p = 0.65; Figure 3C) nor land use type (Kruskal-Wallis; H = 8.79, p = 0.07; Figure 3D).

Table 1

Table 1. Dominant wavelength (λd; mean ± SD) and its trend (Δλd; mean ± SD), summarizing shifts toward shorter λd (i.e., greener or bluer waters) and longer λd (i.e., greener or browner waters) across ecoregions in Missouri.

Figure 4

Figure 4. (A) Annual mean dominant wavelength (λ_d) for all reservoirs (solid line) and grouped by trend category (increasing and decreasing λ_d; dashed lines) from 1984 to 2020. Colors correspond to PDSI index: light blue represents moderately wet years (PDSI > 0), dark blue indicates extremely wet years (PDSI ≥ +3), light brown indicates moderately dry years (PDSI < 0), and dark brown represents extreme drought conditions (PDSI ≤ −3). Graphs (B,D) represent the number of reservoirs distributed across λ_d (binned in 5 nm intervals) in the first and last years between 1984 and 2020; (B) all reservoirs; (C) reservoirs trending toward shorter λ_d; (D) reservoirs trending toward longer λ_d. Horizontal dashed lines in graphs b−d represent the mean λ_d.

Only 6% of reservoirs (n = 10) crossed the 575 nm brown/green threshold, and a single reservoir crossed the 495 nm green/blue threshold. The distribution of all λ_d measurements through time shows that most green reservoirs are closer to the brown rather than blue color threshold (Figure 4B). Across all reservoirs with significant negative (Figure 4C) and positive (Figure 4D) Δλ_d trends, the majority of time series begin and end above 540 nm. Amongst reservoirs trending bluer, the period between the Great Mississippi and Missouri Rivers Floods of 1993 and 2008 constitute the greatest period of change, and more recently λ_d has remained fairly constant, with the exception of wet years when these reservoirs shift back towards browner water (Figure 4A). Amongst reservoirs trending browner, average shifts in water color through time have been comparatively more gradual.

3.2 Drivers of reservoir color change

To quantify potential drivers in water color change (Δλ_d), a Boosted Regression Tree (BRT) analysis was performed with a selection of geomorphological (waterbody elevation, waterbody and watershed area, waterbody and watershed elongation ratio, watershed area to waterbody area ratio) and land cover (percent cover of urban, forest, and agriculture lands in the watershed) predictors, consistent with previous studies (Cao et al., 2023; Liu et al., 2024; Oleksy et al., 2022; Shen et al., 2025). The final BRT model achieved a MAE of 0.24, a RMSE of 0.32, and an R² of 49.9% on independent test data. Overall, waterbody elevation ranked as the most important predictor (RI = 24.7%), followed by urban land cover (RI = 20.3%) and waterbody elongation ratio (RI = 16.3%). The Shapley-derived partial dependence plots for individual predictors (Figure 5) indicate that Missouri reservoirs with increasing λ_d (i.e., browner) are typically located at lower elevations and in areas with a lower percentage of urban land cover (Figures 5A,B). Conversely, reservoirs with decreasing λ_d (i.e., bluer) are found at higher elevations and span a wide range of urban land cover percentages (Figures 5A,B). Morphologically, reservoirs trending browner are generally more elongated, while reservoirs trending bluer are more circular in shape (Figure 5C).

Figure 5

Figure 5. Shapley-based dependence plots for the most important variables selected by the BRT model. Each graph shows how the values of a given predictor (x-axis) influence the model’s prediction of the dominant wavelength slope (Δλ_d), as represented by the SHAP values (y-axis). Positive SHAP values indicate a contribution to increasing Δλ_d (brown dots), while negative values indicate a contribution to decreasing Δλ_d (green dots) across the study reservoirs. Graphs (A−C) correspond to different parameters: (A) waterbody elevation, (B) urban land cover, and (C) waterbody elongation ratio (values close to 1 indicate more circular shapes, while smaller values indicate narrower and more dendritic shapes). RI refers to average relative importance ranked by the selected BRT model. Only significant predictors are shown.

Regardless of overall trends, interannual variations in water color are readily observable across many reservoirs. Averaging across all reservoirs, λ_d varied between 538 and 555 nm through time, where periods of extreme drought and extreme wetness correspond with bluer and browner waters, respectively (Figure 4A). To further explore the interaction between water color and hydroclimate, a linear mixed-effect model including PDSI, ecoregion, their interaction, and year was constructed. Reservoir identity alone (i.e., between-reservoir differences in mean λ_d) accounts for approximately 38% of the explained variance (Table 2). A further 18% of variance is explained by PDSI, ecoregion and their interaction, such that the final model accounts for more than half of the total λ_d variance across the Missouri LimnoSat record (R²c = 0.56; Table 2). PDSI had significant positive effects on λ_d in the Central Irregular Plains, Interior River Valleys and Hills, Ozark Highlands, and Western Corn Belt Plains (Table 3). The strongest positive association was observed in the predominantly agricultural Central Irregular Plains (slope = 1.07, CI = 0.89–1.25). In contrast, the forested, high-relief Ozark Highlands showed a much weaker response in λ_d (slope = 0.19, CI = −0.07–0.45), consistent with reduced erosion during wet periods compared to the low-relief agricultural plains. The only and strongest negative association was observed in the Mississippi Alluvial Plain (slope = −2.97, CI = −9.05–3.12).

Table 2

Table 2. Linear mixed-effect model describing the relationship between dominant wavelength (λ_d; response variable), Palmer Drought Severity Index (PDSI), and ecoregions in Missouri.

Table 3

Table 3. Slopes of Palmer Drought Severity Index (PDSI) on dominant wavelength (λ_d) for each Missouri ecoregion.

3.3 Long-term trends in water quality parameters

The in situ data (minimum of 10 years consecutive data averaged across all reservoirs show increases in Chl a, POM, and the TN:TP ratio over time (Table 4), alongside numerous physical parameters including surface temperature, thermocline depth and buoyancy frequency (i.e., stratification, Table 4). Conversely, PIM and hypolimnetic temperature and dissolved oxygen have significantly decreased over time (Table 4). All other water quality parameters including TP, TN, Secchi depth, and TSS do not have trends. Missouri reservoirs have become warmer and greener (higher Chl a and POM, lower PIM) and more phosphorus-deficient (increased TN:TP ratios), without commensurate trends in water transparency nor absolute nutrient (N and P) concentrations when all in situ measurements are averaged across the state.

Table 4

Table 4. Sen’s slope estimates (mean ± SD) for overall trend and for reservoirs exhibiting increasing or decreasing trends in water quality parameters.

Analysis of reservoir-specific trends; however, provide a more granular description of water quality changes. Of the 130 reservoirs with sufficiently long time series, 108 (83%) had a trend in at least one of the routinely sampled water quality metrics (TP, TN, Chl a, TSS, POM, PIM, Secchi). Single instance trends (i.e., where only a single water quality parameter is changing in a given reservoir) was the most common occurrence (30%), and in these reservoirs, Chl a increases were the most common trend to occur alone (25%), followed by TP (19%). Consistent with statewide trends, Chl a and POM concentrations have increased in 28% and 33% of monitored reservoirs respectively, and decreased in just 5%. PIM has decreased in 24% of reservoirs and increased in 10% (Table 4). TSS, the sum of POM and PIM, has declined (16%) in more reservoirs than where it has increased (11%). Consequently, the water quality indicators with the greatest frequency of significant trends are the POM:TSS and PIM:TSS ratios (Table 4). More than half of all reservoirs have statistically significant trends, with increasing POM:TSS (decreasing PIM:TSS) as the dominant direction (44%, Table 4). With respect to bulk nutrient concentrations, TP has significantly decreased in almost twice as many reservoirs (22%) as it has significantly increased (12%). Moreover, reservoirs with significant decreases in TN and improved water transparency are also more numerous than reservoirs with significant increases in TN and diminished water transparency (Table 4). Thus, at the reservoir scale, more waterbodies are showing improvements in water clarity and declining nutrient concentrations and more reservoirs have increasing Chl a and POM. The most prevalent metric of water quality change is a shift towards a higher fraction of organic matter within the total suspended solids pool (PIM:TSS).

Pairwise linear regression of reservoir-scale water quality slopes, regardless of significance, was conducted on a subset of variables (TP, TN, Chl a, TSS, PIM, POM, Secchi). All slopes were correlated (p < 0.05) and all regression coefficients were positive (e.g., TN and TP slopes positively co-vary), with the exception of water transparency where, unlike the other water quality parameters, a positive slope (deeper transparency) reflects an improvement in water quality. Inspection of the pairwise coefficients of determination (Supplementary Table S8) show that Chl a trends share the strongest association with POM trends (R² = 0.60), followed by TN trends (R² = 0.47). The only other coefficients of determination greater than 0.5 were TN and POM slopes (R² = 0.60), PIM and TSS slopes (R² = 0.64), and POM and TSS slopes (R² = 0.58). Notably, the coefficients of determination between TP and other water quality slopes were small (TN: R² = 0.22, Chl a: R² = 0.12, Secchi: R² = 0.05, Supplementary Table S8).

A second analysis was conducted to quantify how often significant water quality trends co-occur within reservoirs, and separating these into aligned trends (e.g., increasing TP and TN or decreasing TP and TN) and contrasting trends (e.g., increasing TP and decreasing TN). Amongst the reservoirs with increasing Chl a trends, 68% and 0% have increasing and decreasing POM trends, respectively (Supplementary Table S9). For reservoirs with significant TP trends, less than half have commensurate trends in TN, Chl a and other water quality variables, and notably 9% have contrasting TP and Chl a trends (Supplementary Table S9). Overall, with the exception of three pairwise combinations (Chl a and POM, TN and POM, PIM and TSS), all other combinations never exceeded 50%. The pairwise combinations with the greatest frequency of contrasting trends were PIM and Chl a (21%) and POM and PIM (20%).

3.4 Co-variation of water color and water quality parameters

All Missouri water quality variables matched in space and time with LimnoSat data are significantly correlated to λ_d (Spearman rank; p < 0.01; Supplementary Figure S5). Chl a and POM are the least correlated (r = 0.42 and r = 0.42, respectively), and Secchi, TP and PIM are the most correlated (r = −0.67, r = 0.63, and r = 0.63, respectively). To further illustrate these relationships, water quality data were grouped at 10-nm λ_d bins. This analysis yielded statistically significant differences in mean values for all water quality parameters (Kruskal-Wallis; p < 0.05). Subsequent post hoc pairwise analysis was performed across a spectral range (485–575 nm) that encompassed >98% of data and yields one blue group (485–495 nm) and 8 green groups ending at the green-brown transition (Figure 6). Overall, TP increases with λ_d but there are no statistically significantly differences in the 485–555 nm spectral range; however, waterbodies in the 555–565 nm range have significantly higher TP, that in turn are significantly smaller than TP in waterbodies in the 565–575 nm range (Figure 6A). A similar pattern emerges for TN (Figure 6B) and TSS (Figure 6D), with the exception that elevated concentrations of both water quality variables in the two green-brown spectral groups (555–575 nm) are not significantly different from each other. Chl a and POM also increase with λ_d, but unlike all other water quality parameters they reach their respective maxima in the 555–565 nm range then decline at longer wavelengths with no statistically significant differences between groups (Figures 6C–E). PIM follows the same relationship to λ_d as TP, whereby the two brownest waterbody groups are significantly different from each other and from the progressively bluer waterbodies (Figure 6F). Secchi depth is the only water quality parameter to form 4 distinct λ_d groups using the conservative post hoc Bonferroni test (Figure 6G). Similar to TP and PIM, the two brownest spectral end-members constitute two distinct groups, and waterbodies in the 535–555 nm color range are less transparent than bluer (485–535 nm) waterbodies. Finally, DOC has no significant differences across spectral groups but is generally higher in browner, rather than bluer waters (Figure 6H).

Figure 6

Figure 6. Box and whisker plot showing the matchups between in situ water quality parameters and dominant wavelength (λ_d) values binned in 10 nm intervals. Graphs (A−H) correspond to different water quality parameters: (A) total phosphorus (TP), (B) total nitrogen (TN), (C) chlorophyll a (Chl a), (D) total suspended solids (TSS), (E) particulate organic matter (POM), (F) particulate inorganic matter (PIM), (G) Secchi disk depth, and (H) dissolved organic carbon (DOC). Different letters within each graph represent significant differences (p-values ≤ 0.05) among λ_d bins based on pairwise Dunn’s tests with Bonferroni-adjusted p-values. Number in italic below each boxplot represents the number of data available for that bin.

To evaluate the correspondence between water color and water quality trends, we calculated alignment scores (Section 2.3; Table 5) for each reservoir-parameter combination ranging from −1 (consistently opposing significant trends) to 1 (consistently aligned significant trends). Across the 22 reservoirs with statistically significant Δλ_d, all listed water quality variables had positive alignment scores (i.e., alignment by variable in Table 5). PIM and Secchi depth had the highest alignment scores (0.59) and 50% of these reservoirs with significant Δλ_d had statistically significant PIM and Secchi depth trends. POM (0.14) and Chl a (0.07) had the lowest alignment scores. This is broadly consistent with the absence of statistically significant differences across spectral groups, and notable declines in these parameters above 560 nm (Figure 6). Across the 22 reservoirs, 17 had positive alignment scores (i.e., alignment by reservoir in Table 5) such that 77% of reservoirs with significant Δλ_d showed an aggregated change in water quality variables in a manner consistent with the overall correlations to λ_d. However, 5 of the 22 reservoirs (23%) had negative alignment scores that ran counter to expectations. Notably, all negative alignment scores occurred in browning reservoirs (i.e., positive Δλ_d). Amongst these reservoirs, the three with the largest Δλ_d (Hylak_id 1058918, 1057344, 1055671) are all eutrophic systems with mean λ_d values > 540 nm (Supplementary Figures S6–S8). In these reservoirs, Chl a is decreasing, consistent with the observed non-significative decline in this parameter towards the brown end of the spectrum (Figure 6C). Importantly, although these three reservoirs exhibit similar negative alignment scores, they occur in contrasting watershed contexts: Hylak_id 1057344 is located within an urban watershed and has experienced substantial sediment inputs associated with intensive shoreline construction, whereas Hylak_id 1058918 and 1055671 are located in agriculturally dominated watersheds. The fourth reservoir with a negative alignment score (Hylak_id 112285) is also eutrophic and located in an agricultural watershed, though unlike the other reservoirs, PIM shows a modest decline (Supplementary Figure S9). The fifth reservoir (Hylak_id 112606) has the most negative alignment score (−0.86) and is an oligotrophic system located in a forested watershed (Supplementary Figure S10), where small but significant improvements in water quality contrast with a shift toward browner water color.

Table 5

Table 5. Sen’s slope estimates for all reservoirs showing significant changes in water color (Δλ_d) and co-located water quality data.

4 Discussion

The dominant wavelength (λ_d) can be readily derived from a wide array of Earth-observing satellites, and water color is perhaps the most easily understood indicator of water quality among non-experts. In this study, we combine λ_d from the LimnoSat dataset (Topp et al., 2020) with detailed field-based measurements from long-term water quality monitoring programs to investigate changes in both remotely sensed water color and in situ measurements of water quality across Missouri reservoirs. The vast majority of reservoirs across the state are currently green; overall, they are closer to the brown rather than the blue color threshold. Nonetheless, important spatial patterns are present. Reservoirs in the largely forested Ozark Highlands are bluer than those in agriculturally dominated ecoregions, consistent with the statewide spatial gradients in nutrient enrichment previously described (Jones et al., 2008b). Nearly one-third of the study reservoirs displayed trends in water color, of which approximately 60% shifted toward shorter λ_d (n = 91; bluer-green), suggesting that recoveries from more turbid and eutrophic conditions are slightly more common than eutrophication across the LimnoSat record analyzed here (1984–2020). This is consistent with a global scale Δλ_d analysis that showed the majority of inland waterbodies in temperate and boreal regions are also shifting towards shorter wavelengths (Shen et al., 2025).

We further identified climatic, landscape, and waterbody morphological factors driving the spatial and temporal patterns observed. Water color trends in Missouri reservoirs were primarily associated with hydroclimate (PDSI), landscape features (elevation and urban land cover) and reservoir morphology (shape). Similar to other parts of the world, reservoirs across the Midwestern USA are experiencing intensifying hydroclimatic extremes as a result of climate change. For example, climate projections indicate an increased likelihood of drought events in this region, both in frequency and magnitude (Strzepek et al., 2010). Interannual variations in hydrometeorological conditions are significantly correlated with λ_d, such that wet years are associated with higher λ_d values (i.e., browner waters) and drought years are associated with lower λ_d values (i.e., bluer waters). Increased runoff during wet conditions transports greater amounts of terrestrially derived substances, both inorganic (e.g., sediments and nutrients) and organic (e.g., dissolved organic carbon), into reservoirs that collectively shift λ_d to longer wavelengths (De Wit et al., 2016). The effect of hydrological cycles in controlling nutrients (Jones et al., 2008a) and CDOM (Bhattacharya et al., 2022) has been previously explored in Missouri reservoirs, with wet summers increasing watershed-derived inputs. In Missouri, agricultural watersheds typically have greater soil organic matter stores and lower relief compared to forested watersheds (Jones et al., 2008b), conditions that promote higher terrestrial CDOM inputs to reservoirs (Bhattacharya et al., 2022), particularly during periods of high precipitation. This is consistent with continental-scale patterns showing a high prevalence of green and brown eutrophic waterbodies in agricultural Midwest regions (Oliver et al., 2017). Incorporating ecoregions into the linear mixed-effect model improved predictions of the covariation between PDSI and λ_d. Notably, the effect of PDSI on reservoir color is dampened in the forested Ozark Highlands relative to agriculturally dominated ecoregions, consistent with reduced erosion and terrestrial inputs in this high-relief landscape. The only negative association between PDSI and water color was observed in the one reservoir we studied in the Mississippi Alluvial Plain, a pattern likely driven by the extensive hydrologic regulation, levee systems, and floodplain alteration in this ecoregion, which can decouple reservoir optical properties from natural drought/wet cycles. As demonstrated in Missouri reservoirs (Jones et al., 2008a), the absence or presence of water level management across reservoirs is likely an important co-factor modulating the degree to which hydrometeorological variability influences water color.

Reservoir morphology and landscape features collectively explained the spatial variation observed in water color trends. Waterbody elevation consistently ranked as the most important predictor in BRT analysis, revealing that reservoirs trending toward bluer wavelengths are generally situated at higher elevations compared to those trending toward greener wavelengths, consistent with regional analysis (Cao et al., 2023), with the noted exception of high alpine environments (Oleksy et al., 2022). Urban land cover ranked as the second most important predictor in our study. Although urbanization can increase chemical and thermal pollution in freshwater ecosystems (Grimm et al., 2008), reservoirs trended toward shorter λ_d occurred in regions where urban cover ranged from 20% to 80%. Similarly, Oleksy et al. (2022) found that reservoirs shifting towards bluer waters were primarily located in areas with a higher percentage of urban land cover. A potential explanation is that local management practices during summer months may alter water color without necessarily improving overall water quality. In Missouri, for example, several reservoirs are managed to optimize recreational fisheries by herbicide (e.g., glyphosate) applications (Jones et al., 2022). Interventions such as applications of nutrient-containing herbicides can reduce Chl a concentrations by suppressing phytoplankton growth via decreased photosystem II efficiency (Lürling and Roessink, 2006), resulting in short-term shifts toward shorter λ_d. However, over the long-term, these treatments may also increase nutrient availability primarily due to the nitrogen and phosphorus content of these chemicals (Reinl et al., 2022). Reservoir morphology (e.g., waterbody elongation ratio) consistently ranked as the third most important predictor. Reservoirs trending toward shorter λ_d (i.e., bluer waters) are generally smaller and more circular (Figure 5C) and also exhibited greater variability in the magnitude of water color change. This pattern is consistent with findings that smaller waterbodies are more susceptible to shifts in color (Shen et al., 2025).

Changes in water color broadly align with trends in in situ water quality parameters. Co-located λ_d and water quality parameters were significantly related, with reservoirs trending toward shorter λ_d (bluer-green conditions), decreasing TN, TP, PIM, and increasing water transparency being more common than those with opposite trends. Yet, despite this general correspondence, Chl a and POM, both proxies for phytoplankton biomass, increased rather than decreased in more reservoirs. This finding stands in contrast to previous studies that show declines in λ_d are often associated with reductions in Chl a concentrations (Lehmann et al., 2018). Earlier Missouri water quality analyses found that particulate inorganic matter (PIM) in reservoirs is typically the dominant fraction of total suspended solids (TSS), particularly in agricultural watersheds due to erosion and runoff (Jones et al., 2008b; Jones and Knowlton, 2005). However, more recent studies have found that particulate organic matter (POM) concentrations now exceed those of PIM in these same systems (Petty et al., 2020). Our results support this shift, and demonstrate that a shift from PIM to POM is the most commonly occurring metric of water quality change. The majority of significant trends (both negative and positive) are confined within the green spectral region closer to the brown rather than the blue spectral end member, and that in these color regions TP, PIM, and Secchi depth but not Chl a nor POM co-vary with λ_d (Figure 6). Taken together, this demonstrates that diminished λ_d can in fact correspond to elevated phytoplankton biomass, likely through a relaxation of light limitation of phytoplankton communities in eutrophic (i.e., nutrient replete) reservoirs as PIM is much more reflective than POM (Babin et al., 2003). This finding also supports that erosion control resulting from the implementation of best management practices (e.g., shoreline stabilization with rock and water willow; Jones et al., 2022) has resulted in diminished PIM and shifts in water color in three reservoirs.

While significant effort has been undertaken to estimate optically active water quality parameters (e.g., Chl a, TSS, CDOM) from space (Cao et al., 2018; Ondrusek et al., 2012; Smith and Bernard, 2020), this study contributes to a small but growing body of research demonstrating strong linkages between remotely sensed color and non-optically active trophic status indicators, particularly TN and TP (Gao et al., 2015; Guo et al., 2024; Li et al., 2017; Sun et al., 2014; Windle et al., 2025; Xiong et al., 2022). Importantly, these nutrients themselves do not possess optical signatures; instead, their retrieval from satellite data is feasible because of their strong correlation with optically active constituents, which respond to nutrient enrichment through increases in phytoplankton biomass and associated suspended matter (Li et al., 2017). TN and TP concentrations are key chemical indicators of biological growth and eutrophication in freshwater systems (Paerl et al., 2015; 2014); therefore, monitoring these nutrients is essential for effective water quality assessment and management. Future research should further explore advanced modelling techniques, such as machine learning approaches (Zhu et al., 2024), to improve nutrient retrievals from satellite observations and better detect the complex, non-linear relationships between optical properties and nutrient dynamics.

5 Conclusion

In this study, we integrate satellite-derived water color alongside in situ observations, landscape and morphological data to assess long-term water quality trends and identify the factors that drive these trends. Our approach demonstrates the use of remote sensing-derived water color for detecting water quality trends, while underscoring the continued importance of field-based observations to ensure accuracy, validate remote sensing products, and maintain the integrity of water quality assessments. Together, these findings reveal the value of hybrid monitoring frameworks that strategically combine satellite and in situ data.

Shifts in water color reflect changes in watershed land use and reservoir management practices, mediated through changes in particulate inorganic matter, dissolved organic matter, and nutrients. Importantly, our results show that many reservoirs are simultaneously exhibiting decreasing λ_d, reductions in PIM, and deeper Secchi, yet increasing chl a and POM. This finding indicates that historically turbid, agriculturally influenced reservoirs may be transitioning away from strongly light-limited, sediment-dominated conditions toward clearer but still nutrient-rich, phytoplankton-dominated states. Future work should interpret color-based assessments alongside water quality data, particularly nutrient concentrations, especially when informing management decisions, as blueing does not necessarily indicate oligotrophication nor improved trophic conditions. Rather, reductions in suspended mineral particles can relax light limitation, facilitating higher phytoplankton biomass even when water color shifts toward shorter wavelengths.

The observed reduction in brown color reservoirs suggests improvements in water quality in recent years in historically turbid reservoirs. That said, direct water quality data from several of these reservoirs reveals opposing trends in particulate inorganic (decreasing) and organic matter (increasing), underscoring that color changes alone do not necessarily reflect shifts in trophic status. While our study did not directly evaluate the effects of best management practices (BMPs), satellite-derived water color offers a practical tool for assessing how upstream land use and local management practices may affect water quality at the reservoir scale. Previous work in the region has shown that erosion-control BMPs can reduce PIM inputs (Jones et al., 2022), and our findings of long-term declines in PIM and corresponding shifts toward shorter λ_d are consistent with this observation. Future studies that explicitly link BMP implementation histories to long-term water color trends could help quantify management effectiveness at the reservoir scale.

Climate projections indicate increasing frequency and intensity of droughts across the study region (Strzepek et al., 2010). Under such conditions, our findings suggest that reservoirs in the Midwest are likely to shift toward shorter wavelengths, potentially altering primary productivity and energy transfer within these man-made systems. Long-term satellite records of water color therefore offer a promising avenue for anticipating ecosystem responses to climate-driven hydrologic change. In a time of reduced federal funding for long-term water quality monitoring programs, our findings underscore the value of satellite-based assessments as a cost-effective and scalable approach to monitor water quality trends over broad spatial and temporal extents. Many remote-sensing products are now freely available, offering a wide range of spatial (1 m–1 km) and temporal (1–15 days) resolutions, and can therefore provide critical information on water quality trends, even in regions with limited environmental monitoring (Smits et al., 2025). Moreover, monitoring programs can be strategically redesigned to improve data coverage by selecting key site locations and leveraging remote sensing to interpolate across space and time (Smits et al., 2025).

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: The reservoir metadata (IDs and waterbody names) is provided in the Supplementary Material (Supplementary Table S10). The in situ water quality dataset that support the findings of this study are openly available in the Environmental Data Initiative under the following links: https://doi.org/10.6073/pasta/fb18b84f7e5e9d9f23f224e9e6158572 (1984–2018) and https://doi.org/10.6073/pasta/04c5314549c30b280b8dafb9008f8db4 (2019–2020). The water quality profile data are available at https://portal.edirepository.org/nis/mapbrowse?packageid=edi.1544.2.

Author contributions

LP-S: Visualization, Writing – original draft, Data curation, Formal Analysis, Conceptualization, Writing – review and editing, Methodology, Investigation. GS: Investigation, Funding acquisition, Visualization, Writing – review and editing, Formal Analysis, Methodology, Data curation, Supervision, Project administration, Conceptualization. DR: Writing – review and editing, Formal Analysis, Data curation. RN: Funding acquisition, Writing – review and editing, Investigation, Supervision, Data curation, Project administration, Conceptualization.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This research was supported by NASA grant 80NSSC22K0771 (Remote Sensing of Water Quality) awarded to R.L.N and G.M.S. The data is based on work supported by the Missouri Department of Natural Resources (MDNR), which funds the Missouri Statewide Lake Assessment Program (SLAP), coordinated by the University of Missouri (MU) Limnology Laboratory. Part of the data was also obtained through the Lakes of Missouri Volunteer Program (LMVP) funded by MDNR, whose contributions are gratefully acknowledged.

Acknowledgements

The authors are deeply grateful to all students and technicians associated with the MU Limnology Lab over the years, whose efforts in collecting, processing, and analyzing the water quality data made this research possible. We also thank Xiaoxu Guo for reviewing the BRT analysis code and offering valuable feedback on our approach.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2025.1725096/full#supplementary-material

References

Babin, M., Morel, A., Fournier-Sicre, V., Fell, F., and Stramski, D. (2003). Light scattering properties of marine particles in coastal and open ocean waters as related to the particle mass concentration. Limnol. Oceanogr. 48, 843–859. doi:10.4319/lo.2003.48.2.0843

CrossRef Full Text | Google Scholar

Bhattacharya, R., Jones, J. R., Graham, J. L., Obrecht, D. V., Thorpe, A. P., Harlan, J. D., et al. (2022). Nonlinear multidecadal trends in organic matter dynamics in midwest reservoirs are a function of variable hydroclimate. Limnol. Oceanogr. 67, 2531–2546. doi:10.1002/lno.12220

CrossRef Full Text | Google Scholar

Binding, C. E., Jerome, J. H., Bukata, R. P., and Booty, W. G. (2007). Trends in water clarity of the lower Great Lakes from remotely sensed aquatic color. J. Gt. Lakes. Res. 33, 828–841. doi:10.3394/0380-1330(2007)33[828:tiwcot]2.0.co;2

CrossRef Full Text | Google Scholar

Cao, F., Tzortziou, M., Hu, C., Mannino, A., Fichot, C. G., Del Vecchio, R., et al. (2018). Remote sensing retrievals of colored dissolved organic matter and dissolved organic carbon dynamics in North American estuaries and their margins. Remote Sens. Environ. 205, 151–165. doi:10.1016/j.rse.2017.11.014

CrossRef Full Text | Google Scholar

Cao, Z., Melack, J. M., Liu, M., Kutser, T., Duan, H., and Ma, R. (2023). Shifts, trends, and drivers of Lake color across China since the 1980s. Geophys Res. Lett. 50, e2023GL103225. doi:10.1029/2023gl103225

CrossRef Full Text | Google Scholar

Castagna, A., and Vanhellemont, Q. (2025). A generalized physics-based correction for adjacency effects. Appl. Opt. 64, 2719–2743. doi:10.1364/AO.546766

PubMed Abstract | CrossRef Full Text | Google Scholar

Charifson, D. M., Huth, P. C., Thompson, J. E., Angyal, R. K., Flaherty, M. J., and Richardson, D. C. (2015). History of fish presence and absence following Lake acidification and recovery in Lake minnewaska, shawangunk ridge, NY. Northeast Nat. (Steuben) 22, 762–781. doi:10.1656/045.022.0411

CrossRef Full Text | Google Scholar

De Wit, H. A., Valinia, S., Weyhenmeyer, G. A., Futter, M. N., Kortelainen, P., Austnes, K., et al. (2016). Current browning of surface waters will be further promoted by wetter climate. Environ. Sci. Technol. Lett. 3, 430–435. doi:10.1021/acs.estlett.6b00396

CrossRef Full Text | Google Scholar

Dodds, W. K., Bouska, W. W., Eitzmann, J. L., Pilger, T. J., Pitts, K. L., Riley, A. J., et al. (2009). Eutrophication of U. S. Freshwaters: analysis of potential economic damages. Environ. Sci. Technol. 43, 12–19. doi:10.1021/es801217q

PubMed Abstract | CrossRef Full Text | Google Scholar

Elith, J., Leathwick, J. R., and Hastie, T. (2008). A working guide to boosted regression trees. J. Animal Ecol. 77, 802–813. doi:10.1111/j.1365-2656.2008.01390.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Erlandsson, M., Buffam, I., Fölster, J., Laudon, H., Temnerud, J., Weyhenmeyer, G. A., et al. (2008). Thirty-five years of synchrony in the organic matter concentrations of Swedish rivers explained by variation in flow and sulphate. Glob. Chang. Biol. 14, 1191–1198. doi:10.1111/j.1365-2486.2008.01551.x

CrossRef Full Text | Google Scholar

Freeman, C., Evans, C. D., Monteith, D. T., Reynolds, B., and Fenner, N. (2001). Export of organic carbon from peat soils. Nature 412, 785. doi:10.1038/35090628

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Gao, J., Yin, H., Liu, C., Xia, T., Wang, J., et al. (2015). Remote sensing estimation of the total phosphorus concentration in a large lake using band combinations and regional multivariate statistical modeling techniques. J. Environ. Manage 151, 33–43. doi:10.1016/j.jenvman.2014.11.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardner, J. R., Yang, X., Topp, S. N., Ross, M. R. V., Altenau, E. H., and Pavelsky, T. M. (2021). The color of Rivers. Geophys Res. Lett. 48, e2020GL088946. doi:10.1029/2020gl088946

CrossRef Full Text | Google Scholar

Grimm, N. B., Faeth, S. H., Golubiewski, N. E., Redman, C. L., Wu, J., Bai, X., et al. (2008). Global change and the ecology of cities. Sci. (1979) 319, 756–760. doi:10.1126/science.1150195

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, H., Huang, J. J., Zhu, X., Tian, S., and Wang, B. (2024). Spatiotemporal variation reconstruction of total phosphorus in the Great Lakes since 2002 using remote sensing and deep neural network. Water Res. 255, 121493. doi:10.1016/j.watres.2024.121493

PubMed Abstract | CrossRef Full Text | Google Scholar

Hongve, D., Riise, G., and Kristiansen, J. F. (2004). Increased colour and organic acid concentrations in Norwegian forest lakes and drinking water - a result of increased precipitation? Aquat. Sci. 66, 231–238. doi:10.1007/s00027-004-0708-7

CrossRef Full Text | Google Scholar

Jeppesen, E., Søndergaard, M., Jensen, J. P., Havens, K. E., Anneville, O., Carvalho, L., et al. (2005). Lake responses to reduced nutrient loading - an analysis of contemporary long-term data from 35 case studies. Freshw. Biol. 50, 1747–1771. doi:10.1111/j.1365-2427.2005.01415.x

CrossRef Full Text | Google Scholar

Jones, J. R., and Knowlton, M. F. (2005). Suspended solids in Missouri reservoirs in relation to catchment features and internal processes. Water Res. 39, 3629–3635. doi:10.1016/j.watres.2005.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, J. R., Knowlton, M. F., and Obrecht, D. V. (2008a). Role of land cover and hydrology in determining nutrients in mid-continent reservoirs: implications for nutrient criteria and management. Lake Reserv. Manag. 24, 1–9. doi:10.1080/07438140809354045

CrossRef Full Text | Google Scholar

Jones, J. R., Obrecht, D. V., Perkins, B. D., Knowlton, M. F., Thorpe, A. P., Watanabe, S., et al. (2008b). Nutrients, seston, and transparency of Missouri reservoirs and oxbow lakes: an analysis of regional limnology. Lake Reserv. Manag. 24, 155–180. doi:10.1080/07438140809354058

CrossRef Full Text | Google Scholar

Jones, J. R., Obrecht, D., and North, R. L. (2022). Influence of fisheries and shoreline management on limnological characteristics of three Missouri reservoirs. Inland Waters 12, 354–367. doi:10.1080/20442041.2022.2037991

CrossRef Full Text | Google Scholar

King, T. V., Meyer, M. F., Hundt, S. A., Ball, G., Hafen, K. C., Avouris, D. M., et al. (2024). Sentinel-2 ACOLITE-DSF aquatic reflectance for the conterminous United States. U.S. Geol. Surv. Data Release. doi:10.5066/P904243C

CrossRef Full Text | Google Scholar

Leech, D. M., Pollard, A. I., Labou, S. G., and Hampton, S. E. (2018). Fewer blue lakes and more murky lakes across the continental U.S.: implications for planktonic food webs. Limnol. Oceanogr. 63, 2661–2680. doi:10.1002/lno.10967

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehmann, M. K., Nguyen, U., Allan, M., and van der Woerd, H. J. (2018). Colour classification of 1486 lakes across a wide range of optical water types. Remote Sens. (Basel) 10, 1273. doi:10.3390/rs10081273

CrossRef Full Text | Google Scholar

Li, Y., Zhang, Y., Shi, K., Zhu, G., Zhou, Y., Zhang, Y., et al. (2017). Monitoring spatiotemporal variations in nutrients in a large drinking water reservoir and their relationships with hydrological and meteorological conditions based on landsat 8 imagery. Sci. Total Environ. 599–600, 1705–1717. doi:10.1016/j.scitotenv.2017.05.075

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Kim, S., Glamore, W., Tamburic, B., and Johnson, F. (2024). Remote sensing of water colour in small southeastern Australian waterbodies. J. Environ. Manage 352, 120096. doi:10.1016/j.jenvman.2024.120096

PubMed Abstract | CrossRef Full Text | Google Scholar

Lürling, M., and Roessink, I. (2006). On the way to cyanobacterial blooms: impact of the herbicide metribuzin on the competition between a green alga (scenedesmus) and a cyanobacterium (microcystis). Chemosphere 65, 618–626. doi:10.1016/j.chemosphere.2006.01.073

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, M. F., Topp, S. N., King, T. V., Ladwig, R., Pilla, R. M., Dugan, H. A., et al. (2024). National-scale remotely sensed lake trophic state from 1984 through 2020. Sci. Data 11, 77. doi:10.1038/s41597-024-02921-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Monteith, D. T., Stoddard, J. L., Evans, C. D., De Wit, H. A., Forsius, M., Høgåsen, T., et al. (2007). Dissolved organic carbon trends resulting from changes in atmospheric deposition chemistry. Nature 450, 537–540. doi:10.1038/nature06316

PubMed Abstract | CrossRef Full Text | Google Scholar

Mouw, C. B., Greb, S., Aurin, D., DiGiacomo, P. M., Lee, Z., Twardowski, M., et al. (2015). Aquatic color radiometry remote sensing of coastal and inland waters: challenges and recommendations for future satellite missions. Remote Sens. Environ. 160, 15–30. doi:10.1016/j.rse.2015.02.001

CrossRef Full Text | Google Scholar

Nakagawa, S., and Schielzeth, H. (2013). A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods Ecol. Evol. 4, 133–142. doi:10.1111/j.2041-210x.2012.00261.x

CrossRef Full Text | Google Scholar

Oleksy, I. A., Collins, S. M., Sillen, S. J., Topp, S. N., Austin, M., Hall, E. K., et al. (2022). Heterogenous controls on lake color and trends across the high-elevation U.S. rocky Mountain region. Environ. Res. Lett. 17, 104041. doi:10.1088/1748-9326/ac939c

CrossRef Full Text | Google Scholar

Oleksy, I. A., Solomon, C. T., Jones, S. E., Olson, C., Bertolet, B. L., Adrian, R., et al. (2024). Controls on Lake pelagic primary productivity: formalizing the nutrient-color paradigm. J. Geophys Res. Biogeosci 129, e2024JG008140. doi:10.1029/2024jg008140

CrossRef Full Text | Google Scholar

Oliver, S. K., Collins, S. M., Soranno, P. A., Wagner, T., Stanley, E. H., Jones, J. R., et al. (2017). Unexpected stasis in a changing world: lake nutrient and chlorophyll trends since 1990. Glob. Chang. Biol. 23, 5455–5467. doi:10.1111/gcb.13810

PubMed Abstract | CrossRef Full Text | Google Scholar

Ondrusek, M., Stengel, E., Kinkade, C. S., Vogel, R. L., Keegstra, P., Hunter, C., et al. (2012). The development of a new optical total suspended matter algorithm for the Chesapeake Bay. Remote Sens. Environ. 119, 243–254. doi:10.1016/j.rse.2011.12.018

CrossRef Full Text | Google Scholar

Paerl, H. W., Xu, H., Hall, N. S., Zhu, G., Qin, B., Wu, Y., et al. (2014). Controlling cyanobacterial blooms in hypertrophic Lake Taihu, China: will nitrogen reductions cause replacement of non-N2 fixing by N2 fixing taxa? PLoS One 9, e113123. doi:10.1371/journal.pone.0113123

PubMed Abstract | CrossRef Full Text | Google Scholar

Paerl, H. W., Xu, H., Hall, N. S., Rossignol, K. L., Joyner, A. R., Zhu, G., et al. (2015). Nutrient limitation dynamics examined on a multi-annual scale in Lake Taihu, China: implications for controlling eutrophication and harmful algal blooms. J. Freshw. Ecol. 30, 5–24. doi:10.1080/02705060.2014.994047

CrossRef Full Text | Google Scholar

Petty, E. L., Obrecht, D. V., and North, R. L. (2020). Filling in the flyover zone: high phosphorus in midwestern (USA) reservoirs results in high phytoplankton biomass but not high primary productivity. Front. Environ. Sci. 8, 111. doi:10.3389/fenvs.2020.00111

CrossRef Full Text | Google Scholar

R Core Team (2024). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Available online at: https://www.R-project.org/.

Google Scholar

Reinl, K. L., Harris, T. D., Elfferich, I., Coker, A., Zhan, Q., De Senerpont Domis, L. N., et al. (2022). The role of organic nutrients in structuring freshwater phytoplankton communities in a rapidly changing world. Water Res. 219, 118573. doi:10.1016/j.watres.2022.118573

PubMed Abstract | CrossRef Full Text | Google Scholar

Roulet, N., and Moore, T. R. (2006). Browning the waters. Nature 444, 283–284. doi:10.1038/444283a

PubMed Abstract | CrossRef Full Text | Google Scholar

Schindler, D. W., Hecky, R. E., and McCullough, G. K. (2012). The rapid eutrophication of Lake Winnipeg: greening under global change. J. Gt. Lakes. Res. 38, 6–13. doi:10.1016/j.jglr.2012.04.003

CrossRef Full Text | Google Scholar

Shen, X., Ke, C. Q., Duan, Z., Cai, Y., Li, H., and Xiao, Y. (2025). Satellite observations reveal widespread color variations in global Lakes since the 1980s. Water Resour. Res. 61, e2023WR036926. doi:10.1029/2023wr036926

CrossRef Full Text | Google Scholar

Sillen, S. J., Ross, M. R. V., and Collins, S. M. (2024). Long-term trends in productivity across intermountain west Lakes provide no evidence of widespread eutrophication. Water Resour. Res. 60, e2023WR034997. doi:10.1029/2023wr034997

CrossRef Full Text | Google Scholar

Smith, M. E., and Bernard, S. (2020). Satellite ocean color based harmful algal bloom indicators for aquaculture decision support in the southern Benguela. Front. Mar. Sci. 7, 61. doi:10.3389/fmars.2020.00061

CrossRef Full Text | Google Scholar

Smith, B., Pahlevan, N., Schalles, J., Ruberg, S., Errera, R., Ma, R., et al. (2020). A Chlorophyll-a algorithm for Landsat-8 based on mixture density networks. Front. Remote Sens. 1, 623678. doi:10.3389/frsen.2020.623678

CrossRef Full Text | Google Scholar

Smits, A. P., Hall, E. K., Deemer, B. R., Scordo, F., Barbosa, C. C., Carlson, S. M., et al. (2025). Too much and not enough data: challenges and solutions for generating information in freshwater research and monitoring. Ecosphere 16, e70205. doi:10.1002/ecs2.70205

CrossRef Full Text | Google Scholar

Strzepek, K., Yohe, G., Neumann, J., and Boehlert, B. (2010). Characterizing changes in drought risk for the United States from climate change. Environ. Res. Lett. 5, 044012. doi:10.1088/1748-9326/5/4/044012

CrossRef Full Text | Google Scholar

Sukristiyanti, S., Maria, R., and Lestiana, H. (2018). “Watershed-based morphometric analysis: a review,” in IOP conference series: earth and environmental science (Bristol, England: IOP Publishing).

Google Scholar

Sun, D., Qiu, Z., Li, Y., Shi, K., and Gong, S. (2014). Detection of total phosphorus concentrations of turbid inland waters using a remote sensing method. Water Air Soil Pollut. 225, 1953. doi:10.1007/s11270-014-1953-6

CrossRef Full Text | Google Scholar

Topp, S., Pavelsky, T., Yang, X., Gardner, J., and Ross, M. R. V. (2020). LimnoSat-US: a remote sensing dataset for U.S. Lakes from 1984-2020. Zenodo. doi:10.5281/zenodo.4139695

CrossRef Full Text | Google Scholar

Topp, S. N., Pavelsky, T. M., Dugan, H. A., Yang, X., Gardner, J., and Ross, M. R. V. (2021). Shifting patterns of summer Lake color phenology in over 26,000 US Lakes. Water Resour. Res. 57, e2020WR029123. doi:10.1029/2020WR029123

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Li, J., Shen, Q., Zhang, B., Zhang, F., Lu, Z., et al. (2014). “MODIS-based radiometric color extraction and classification of inland water with the Forel-Ule scale: A case study of Lake Taihu,” in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. 8 (2), 907–918. doi:10.1109/JSTARS.2014.2360564

CrossRef Full Text | Google Scholar

Webster, K. E., Soranno, P. A., Cheruvelil, K. S., Bremigan, M. T., Downing, J. A., Vaux, P. D., et al. (2008). An empirical evaluation of the nutrient-color paradigm for lakes. Limnol. Oceanogr. 53, 1137–1148. doi:10.4319/lo.2008.53.3.1137

CrossRef Full Text | Google Scholar

Williamson, C. E., Morris, D. P., Pace, M. L., and Olson, O. G. (1999). Dissolved organic carbon and nutrients as regulators of lake ecosystems: resurrection of a more integrated paradigm. Limnol. Oceanogr. 44, 795–803. doi:10.4319/lo.1999.44.3_part_2.0795

CrossRef Full Text | Google Scholar

Windle, A. E., Malkin, S. Y., Hood, R. R., and Silsbe, G. M. (2025). Optical water typing in optically complex waters: a case study of Chesapeake Bay. Sci. Total Environ. 981, 179558. doi:10.1016/j.scitotenv.2025.179558

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, J., Lin, C., Cao, Z., Hu, M., Xue, K., Chen, X., et al. (2022). Development of remote sensing algorithm for total phosphorus concentration in eutrophic lakes: conventional or machine learning? Water Res. 215, 118213. doi:10.1016/j.watres.2022.118213

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., O’Reilly, C. M., Gardner, J. R., Ross, M. R. V., Topp, S. N., Wang, J., et al. (2022). The color of earth’s Lakes. Geophys Res. Lett. 49, e2022GL098925. doi:10.1029/2022gl098925

CrossRef Full Text | Google Scholar

Zhu, L., Cui, T., Runa, A., Pan, X., Zhao, W., Xiang, J., et al. (2024). Robust remote sensing retrieval of key eutrophication indicators in coastal waters based on explainable machine learning. ISPRS J. Photogrammetry Remote Sens. 211, 262–280. doi:10.1016/j.isprsjprs.2024.04.007

CrossRef Full Text | Google Scholar