Data used in deriving bbp from CALIPSO comes from NASA/CNES’s LiDAR Level 1B profile data, Version 4–10 product (Winker et al., 2009). For the majority of this assessment, we followed Behrenfeld et al. (2013), implementing changes that have come about over the last decade to improve the reliability of the results (Lu et al., 2013; 2014; 2021b; Bisson et al., 2021). A schematic of how b _bp is derived from CALIPSO is shown in Figure 1. For the scope of this analysis, b _bp refers to the backscatter sampled at 532 nm. At every point along the satellite track, the co-polarized and cross-polarized channel returns are extracted. A cross-talk correction between the two channels is implemented and the transient response from the surface is removed (Lu et al., 2014; 2021b). The corrected signal is then used to calculate a depolarization ratio (δ _t ) between the two channels for the first three bins below the surface of the water. Following this, a series of filtering is done to eliminate signals that would result in a contaminated result. Implementation of this filtering is as follows, as was done in Dionisi et al. (2020):

1) Removal of signal that is, flagged as saturated in the data product.

Lu et al. (2018) implemented a signal saturation flag to the CALIPSO data product. Here, we only consider data that is, not saturated in any way, and ignore data that is, flagged as possibly saturated or certainly saturated as this signal would not yield reliable results.

2) Removal of signal that had cloud coverage.

Clouds are identified though two sources. First, if the water surface peak from the LiDAR return is not within 120 m of the actual surface (derived from the Digital Elevation Model (DEM) flag on the CALIPSO data), then the signal is considered to be polluted by clouds. Second, if the integrated attenuated backscatter (IAB) for the entire LiDAR return is greater than a threshold value (0.017sr ⁻1), then the signal is considered contaminated by clouds (Dionisi et al., 2020).

3) Removal of the signal where the depolarization ratio (δ _t ) exceeded 0.5.

Realistic values of the depolarization ratio (δ _t ) certainly would be below 0.5 (Dionisi et al., 2020). As such, all data returns with a depolarization ratio greater than 0.5 are not considered.

4) Removal of the signal where the wind speed is less than 2 m/s and greater than 9 m/s.

Low wind speeds (less than 2 m/s) result in signal saturation and high wind speeds (greater than 9 m/s) result in turbid waters. As each of these cases would result in an unreliable signal, they are therefore not considered.

After preliminary filtering of the signal, we were able to start the calculation of b _bp . This was done through the use of parameters taken from Behrenfeld et al. (2013), Bisson et al. (2021), and though two dynamic variables. A listing of the constants and values are shown in Table 1. The first dynamic variable used in deriving b _bp is the diffuse attenuation coefficient for downwelling irradiance (K _d ), which is obtained from MODIS optical imagery. Specifics surrounding the acquisition of K _d are shown in Section 2.2. We chose to directly use the MODIS derived K _d measurements rather than using the empirical relationship for K _d in Bisson et al. (2021) because we are analyzing a freshwater environment. As such, the relationship between MODIS K _d and the depolarization ratio (δ _t ) may be different. However, it is likely that using either method will result in a very similar result, as the K _d used in Bisson et al. (2021) is still derived from MODIS via an empirical relationship.

The second and most important dynamic variable in the derivation of b _bp from CALIPSO is water surface wind speed (v). Wind speed is used in deriving wave height (Cox and Munk, 1954; Hu et al., 2008), which is directly used in calculating b _bp from the depolarization ratio (δ _t ). For every point along the CALIPSO tracks, we used a dynamic wind speed obtained from NOAA GLCFS (NOAA, 2022). The high temporal resolution of the wind speed model allowed us to have wind speed measurements down to the same hour of each CALIPSO flyover. High winds will result in waves that make the water too turbid to obtain reliable b _bp measurements and low wind speeds can result in signal saturation (Behrenfeld et al., 2013). As such, we implemented a filter by removing all measurements that had a wind speed greater than 9 m/s and less than 2 m/s (Dionisi et al., 2020), as is shown in the pre-processing steps. Thus, we can now define b _bp as a function of the depolarization ratio (δ _t ), wind speed v), the attenuation coefficient for downwelling irradiance (K _d ), and the combination of previously defined constants C) following the relationship shown in Behrenfeld et al. (2013) as: b b p = C ∗ f δ t , v , K d (1)

Level 2 MODIS imagery intersecting Lake Michigan was downloaded through the NASA Ocean Biology Processing Group (OBPG; https://oceancolor.gsfc.nasa.gov/). Each image was processed using the Color Producing Agents Algorithm (CPA-A; Shuchman et al. (2013)) in order to derive estimates of chlorophyll-a concentration, suspended minerals concentration, and CDOM (Colored Dissolved Organic Matter) absorption. Using these three estimates, bulk absorption and bulk backscatter coefficients were derived for the following MODIS bands: 412, 443, 488, 531, 547, and 667 nm. At the same bands, b _bp was computed from the bulk backscatter by removing the backscatter due to pure water (coefficients derived from Morel et al. (1974). The diffuse attenuation coefficient (K _d ) at the above wavelengths was estimated using a method outlined in Lee et al. (2005), which uses the bulk absorption and backscatter coefficients as well as the solar zenith angle.

Yearly average images were also computed for both b _bp and K _d . First, daily average images were generated by computing the mean of overlapping pixels within all satellite images from a given day. The yearly average images were then computed as the mean of all daily images within that year.

In situ b _bp was sampled by the National Oceanic and Atmospheric Administration’s Great Lakes Environmental Research Laboratory (NOAA GLERL). This was done primarily in the spring (March-May) and summer (June-August), with a scattering of samples in the fall (September- November), at several stations on Lake Michigan between 2015 and 2019. Observations of b _bp were derived from data collected by a WET Labs BB9 sensor, which measures volume scattering coefficients at 9 wavelengths (412, 440, 488, 510, 532, 595, 650, 676, and 715 nm). During sampling, the BB9 is mounted in a package along with other sensors including a WET Labs ac-s, Sea-Bird CTD, and WET Labs fluorometer. Packaging these sensors provides concurrent measurements of salinity, temperature, and absorption which are necessary for processing BB9 data to b _bp . The package was deployed vertically through the water column using a crane.

Using the WAP software package (WET Labs), ac-s, CTD, and BB9 data were converted from binary data to text files, and BB9 data were processed to b _bp using protocols outlined in Zaneveld et al. (2003). First, the total volume scattering function (β _t ) is corrected using the coincident total absorption (a _t ) measurements from the ac-s after having been re-sampled to the BB9 wavelengths. Next, the volume scattering function of the water (β _w ) was calculated according to Boss and Pegau (2001), utilizing the coincident CTD-measured temperature and salinity. The particulate fraction of the volume scattering function (β _p ) is calculated as the difference between β _t and β _w . b _bp is then computed according to the following equation using a χ factor of 1.1 (Sullivan et al., 2013): b b p = 2 π ∗ χ ∗ β p (2) Finally, the b _bp data is binned to 1 m with the vertical profiles then averaged between 0 and 50 m below the water surface.

For the scope of this assessment, we chose to limit our study to only Lake Michigan rather than any of the other Great Lakes. This was done purposefully for a two main reasons. First, the way the CALISPO flyovers were oriented coincided very well with the geometry of the lake. For Lake Michigan, the satellite had two unbroken and intersecting day/night tracks that spanned a few degrees of latitude (Figure 2). This match up allowed us to effectively preform our analysis even with the limited spatial coverage of the CALISPO satellite. Secondly, the distribution of in situ sampled b _bp values was the highest in Lake Michigan. This distribution of samples along similar lines of latitude to that of the CALIPSO tracks allowed us to compare our derived product effectively. A map of both of the CALISPO tracks used in this survey along with the locations of all in situ sampling stations is shown in Figure 2.

As a first test of the ability to derive b _bp from CALIPSO, we computed an average b _bp across the entire lake for every year in the data record. We did this for both the daytime and nighttime CALIPSO tracks. To ascertain how this value compared to other measurements of b _bp , we took a yearly average b _bp for the MODIS data at the same location as the CALIPSO data. In conjunction with CALIPSO and MODIS, we also calculated a yearly average b _bp for the in situ data across the entire lake. We compared the three metrics in Figure 3A (daytime) and Figure 3B (nighttime). As MODIS is unable to sample at nighttime and as there is no documented nighttime b _bp samples on Lake Michigan, the nighttime b _bp from CALIPSO is compared to daytime measurements.

Our results for the daytime flyovers showed more yearly variability in the CALIPSO b _bp than the MODIS b _bp (Figure 3A). However, over the course of the 15 year period, there was no discernible trend in b _bp (p-value > 0.05 ). The 95% confidence bounds for the CALIPSO b _bp were similar to that of the MODIS derived results, and much smaller than the bounds on the in situ sampling. These results also held for the nighttime CALIPSO results (Figure 3B). However, the nighttime results are systematically lower for every year in the record when compared to their daytime counterparts and the MODIS/in situ data.

Both cloud cover and high wind speeds limited the return rate of usable CALIPSO data and therefore did not allow us resolve seasonal trends across the lake on a yearly basis. However, because trends across the 15 years time period of the CALIPSO, MODIS, and in situ data were largely unchanged (relative to the standard error of the in situ measurements), we felt justified in combining the entire time record into a seasonally divided data set and then evaluating this data set spatially across the lake. We did this for both the daytime and nighttime measurements. To start, the daytime measurements across Lake Michigan are shown in Figure 4A. These results are divided up into four seasons: Spring (March through May), Summer (June though August), Fall (September though November) and Winter (December through February).

At a first order evaluation, for the spring, summer, and fall we see very good coherence between all three methods of collecting Daytime b _bp (Figure 4A). In the winter, we see a much larger divergence between measurements, which is not surprising in that the MODIS result is not well calibrated for winter. This could be related to one or more of the following potential issues. First, the CPA-A, which is used to generate the MODIS b _bp estimates, is calibrated based on in situ measurements throughout the lakes (Shuchman et al., 2013). However, because our in situ dataset does not include any measurements collected in winter, it is unclear how suitable the current calibration is during the winter season. Second, the CALIPSO result may also be influenced by ice in the lake. With these issues in mind, CALIPSO may be able to supplement non-existent winter sampling campaigns if validation of the result could be performed. The nighttime measurements from CALILPSO again show systematically lower response across all sections of the lake when compared to the daytime sampled results (Figure 4B). The usable data retrieval rate of the nighttime measurements was also nearly an order of magnitude higher then that of the daytime (7% vs. 1%).

Also of note is the lack of CALIPSO daytime data between 42.5° and 43.5°latitude. This likely is a direct result of the optical complexity of the waters of Lake Michigan in this region. Satellite optical imagery frequently shows the existence of sediment plumes in this part of the lake (Lohrenz et al., 2004; Vanderploeg et al., 2007), which may be resulting in a CALIPSO return that is, flagged as contaminated (for one or more of the previously shown filtering steps). Likewise, this part of the CALIPSO track is mostly nearshore, which further increases the optical complexity of the water and may result in further signal loss. This is also related to the large variability of in situ values in these optically complex waters, which are likely to have more variability in their b _bp relative to portions of the lake that are more spatially consistent.

To a higher level analysis of the daytime results, there is some smaller scale divergence across the lake between the CALIPSO and MODIS measurements. This is especially prevalent in the spring time (Figure 4A). Taking the in situ values as being ground truth, we next compared the daytime CALIPSO and MODIS b _bp results to the their closest measurement spatially on a seasonal basis, taking a median percent error for each season and instrument. We did this for both the spring and the summer separately (when in situ measurements were available) and then also combined the results across all seasons (Figure 5).

We found that CALIPSO derived b _bp showed better agreement relative to the in situ sampling than the MODIS derived b _bp . This was especially true in the springtime where the CALIPSO b _bp had a median percent error of 7% and the MODIS b _bp had a median percent error of 27%. In the summer, both instruments were nearly the same in their performance, with the CALIPSO (14% error) only sightly outperforming the MODIS (15% error). Finally, taken together regardless of season, CALIPSO (8%) was closer to in situ values than MODIS (18%).

Our results indicate that CALIPSO can retrieve a reliable b _bp signal on a large, freshwater lake. However, deriving trends with a higher resolution than yearly across the entire lake or seasonally across the entire data record is impractical using the CALIPSO data. This is due to a weather dependent return rate of usable CALIPSO backscatter data across the lake. For daytime measurements, the amount of good measurements after filtering is around 1%. This improves substantially for the nighttime measurements where the amount of usable data climbs to approximately 7%. However, even at 7% retrieval, the limited spatially coverage of CALIPSO prevents a more in depth analysis.

Reasons for the low usable data percentage of CALIPSO on the Great Lakes stems mostly from two sources. First, average wind speeds on Lake Michigan are generally around 6 m/s, with values varying both spatially and temporally (Li et al., 2010). Due to turbid waters at high wind speeds and to signal saturation at low wind speeds, a range of wind speeds of between 2 m/s and 9 m/s is required in order to reliable derive b _bp from CALIPSO (Behrenfeld et al., 2013). Many of the CALIPSO flyovers on Lake Michigan take place when the wind speed is greater than the maximum allowed wind speed, resulting in a considerable loss of data. The second major source of data loss comes from the cloud coverage on the great lakes, where the percentage of cloud free days is less than 50% (Ju and Roy, 2008). Clouds prevent reliable retrieval of the signal from CALIPSO and therefor result in a null measurement.

A final note on the return rate of usable data from CALIPSO is the substantially higher retrieval rate in the nighttime hours to that of the daytime hours (7% vs. 1%). This is likely due to the behavior of clouds on wind speeds on Lake Michigan between the daytime and the nighttime. In the daytime, temperature gradients between the lake and the land produce high winds, an effect which may be diminished in the nighttime when temperature gradients are much less steep (Laird et al., 2001). This could result in both less cloud coverage and lower wind speeds on the lake in the evening hours, resulting in a higher data usability rate.

As the usable data retrieval rate for nighttime measurements is higher for CALIPSO, it would be advantageous to use nighttime measurements of b _bp to further monitor the Great Lakes. However, our results indicate that there is a substantial offset in nighttime b _bp across all years (Figure 3B) and seasons (Figure 4B). This offset is sometimes more than 50% lower than the closest daytime measurement. Theory on the open ocean suggests that the nighttime measurement should be intrinsically 10% lower due to the diurnal size differences in particulates (Kheireddine and Antoine, 2014; Behrenfeld et al., 2019). This difference could be exacerbated in the freshwater ecosystem where zooplankton and phytoplankton are stoichiometricly distinct compared to their marine equivalent, and therefore could have much different and more amplified diurnal difference to their combined back-scattering (Elser and Hassett, 1994).

The challenge with validating the nighttime CALIPSO measurements is the lack of other sources of nighttime b _bp to compare it to. MODIS can only sample in the daytime as it is an optical instrument and there has not yet been any effort to sample b _bp in the nighttime across the Great Lakes. Because all current and future LiDAR based satellites will sample in the nighttime, and because the retrieval rate of usable data for nighttime measurements is nearly an order of magnitude better than daytime measurements, future sampling efforts across the Great Lakes should be gauged to have a nighttime component. This would serve to validate any future spaceborne LiDAR derived b _bp measurements.

Our results indicated that daytime derived CALIPSO b _bp aligned better with the in situ sampling when compared to MODIS derived b _bp . This was especially true in the springtime where CALIPSO measurements were more than 20% closer to the in situ values than MODIS measurements. However, in the summer months CALIPSO was only slightly closer than MODIS, with a difference between them of less than 1%. This result is in line with previous studies on the global oceans, where CALIPSO performed better than MODIS when compared to in situ data gathered by the network of ARGO floats (Bisson et al., 2021). However, it should be noted that in situ sampling is quite variable and further analysis would be needed to further examine the performance of the CALIPSO measurements to MODIS measurements. For example, in situ measurements of b _bp on Lake Michigan are taken only periodically (usually twice a year) and only at one particular section of the lake. These sampling campaigns also are done via shipborne collection, which are intrinsically time consuming. To set up a more consistent and more efficient sampling campaign, it would be advantageous to establish a system of floats (similar to the open ocean) that could collect in situ values of b _bp regularly throughout the year. This would vastly improve the analysis.

The difference in reliability between the spring and summer for MODIS is likely due the summer biasing of the b _bp derivation from optical imagery. MODIS derivered b _bp is calculated, in part, by using in situ values to calibrate the method. Most of the in situ sampling that is, used to calibrate the MODIS derived product comes from summertime measurements. This results in a heavy biasing towards the summer months which yields a summertime MODIS b _bp that aligns better with in situ measurements and a springtime MODIS b _bp that diverges. Moreover, there is very little difference in coherence between seasons for CALIPSO because CALISPO derived b _bp is independent of in situ sampling campaigns.

Here, we derived b _bp from a spaceborne, LiDAR based, satellite on a large freshwater lake. We found that b _bp derived in this manner matches well with in situ sampled and MODIS derived b _bp values. We also found that the LiDAR derived values tend to be closer than the MODIS derived values when compared to the in situ values, however variability in the in situ sampling may be biasing this relationship. That said, the practicality of CALIPSO derived b _bp is limited on the Great Lakes due to three main reasons:

1) The weather dependent retrieval rate of daytime measurements is less than 1%, which makes monitoring small scale trends nearly impossible.

2) The spatial coverage of CALIPSO is limited in the scope of the Great Lakes, where the satellite only makes a few flyovers across repeat tracks.

3) CALIPSO, after 15 years in service, is nearing the end of its usable life and therefor further data that will be acquired by the satellite is likely minimal.

With these drawbacks in mind, the usability of CALIPSO derived b _bp on the Great Lakes likely lies in it ability to supplement in situ measurements, which are used to validate the gridded b _bp MODIS products. Previously, we stated that CALISPO was more in line with in situ measurements than MODIS, especially in the springtime. This is due to summer biasing which is related to the heavy summer distribution of in it situ measurements. However, CALIPSO derived results may be able to serve as proxy “in situ” b _bp values. This would greatly supplement current sampling efforts and improve MODIS derived b _bp products.

Even with CALISPO coming to an end, the future of spaceborne LiDAR derived b _bp on the Great Lakes is still bright. Recent studies on the global have used a new LiDAR based satellite that was launched in 2018, ICESat-2, to calculate b _bp both as an along track variable and as a function of water depth (Lu et al., 2020; 2021b; a). With considerably higher spatial coverage than CALIPSO and the ability to profile b _bp at depth, ICESat-2 could provide valuable information about water quality on the Great Lakes. With that in mind, we believe that spaceborne LiDAR will be a major component of monitoring efforts on the Great Lakes over the next 10 years.