Carbon Fixation Trends in Eleven of the World’s Largest Lakes: 2003–2018

: Large freshwater lakes provide immense value to the surrounding populations, yet there is limited understanding of how these lakes will respond to climate change and other factors. This study uses satellite remote sensing to estimate annual, lake-wide primary production in 11 of the world’s largest lakes from 2003–2018. These lakes include the ﬁve Laurentian Great Lakes, the three African Great Lakes, Lake Baikal, and Great Bear and Great Slave Lakes. Mean annual production in these lakes ranged from under 200 mgC / m 2 / day to over 1100 mgC / m 2 / day, and the lakes were placed into one of three distinct groups (oligotrophic, mesotrophic, or eutrophic) based on their level of production. The analysis revealed only three lakes with signiﬁcant production trends over the study period, with increases in Great Bear Lake (24% increase over the study period) and Great Slave Lake (27%) and a decline in Lake Tanganyika ( − 16%). These changes appear to be related to climate change, including increasing temperatures and solar radiation and decreasing wind speeds. This study is the ﬁrst to use consistent methodology to study primary production in the world’s largest lakes, allowing for these novel between-lake comparisons and assessment of inter-annual trends. remote sensing-derived monthly mean data were used to estimate the light attenuation coe ﬃ cient in the visible domain (350–700 nm), K PAR , and photic zone depth (the depth where 1% of surface light remains), z 1% , using a methodology outlined in Lee This approach models light attenuation as a function of the absorption and backscattering coe ﬃ cients of the water (a, and b b ) and solar zenith angle. Hourly values of K PAR and z 1% were derived using the monthly averaged coe ﬃ cients and hourly solar zenith angle was determined based on the Julian date and latitude [49]. Latitude was set to be the approximate geometric center of each lake.


Introduction
Freshwater lakes make up less than 1% of the total global water volume [1], yet they provide tremendous value to the planet and the nearby human populations. While the number of global lakes has been estimated to be in the hundreds of millions [2,3], lakes are unevenly distributed in size, with the five largest lakes containing more than half the planet's liquid fresh water [4]. These large lakes provide many important services to the surrounding human populations, both socially and economically. This is particularly true for the lakes in more populous regions, including the African Great Lakes (Lakes Malawi, Tanganyika, and Victoria) and Laurentian Great Lakes (Lakes Erie, Huron, Michigan, Ontario, and Superior). For instance, more than 4.5 billion m 3 of water is extracted annually from the Laurentian Great Lakes for public use [4], including drinking water for approximately 35 million people [5]. Commercial fishing on the large lakes provides both sustenance and employment to nearby populations. The tropical African Great Lakes have the largest fisheries, with a total annual catch of nearly 1.3 megatons [4], providing a bulk of the nutrition for the local residents [4,6]. Large lakes also provide tremendous economic value through recreation and tourism; this industry has been estimated to bring in over USD 6 billion dollars annually in the Laurentian Great Lakes [4]. Other services provided by the large lakes include transportation, hydropower, and irrigation [4][5][6]. In addition to the clear social and economic value these lakes provide to the local populations, these lakes are also an important part of the global carbon cycle [7][8][9], however their total contribution in terms of carbon fixation is not well understood.
will allow for the first direct comparisons between these large lakes as well as time series analyses to assess the impacts of recent climate change on lake-wide primary production.

Materials and Methods
Primary production was estimated for 11 of the world's largest lakes ( Figure 1) using remotely sensed chlorophyll-a concentration, water temperature, and light extinction coefficients; empirical relationships between water temperature and P max (maximum photosynthetic rate at light saturation) and modeled surface irradiance. This approach is similar to that used by Fahnenstiel et al. [38] with several modifications.
Water 2020, 12, x FOR PEER REVIEW 3 of 16

Materials and Methods
Primary production was estimated for 11 of the world's largest lakes ( Figure 1) using remotely sensed chlorophyll-a concentration, water temperature, and light extinction coefficients; empirical relationships between water temperature and Pmax (maximum photosynthetic rate at light saturation) and modeled surface irradiance. This approach is similar to that used by Fahnenstiel et al. [38] with several modifications. The Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on board the Aqua satellite was the source of remotely sensed imagery. This imagery was acquired through the National Aeronautics and Space Administration (NASA) Ocean Biology Processing Group (OBPG) (https://oceancolor.gsfc.nasa.gov/) at Level 2, having been processed by the OBPG using their standard atmospheric correction routines [46]. For each lake, all images collected between 1 January 2003 and 31 December 2018 that intersected any portion of the target were identified and downloaded. Because MODIS Aqua is a polar-orbiting sensor, lakes closer to the poles are observed more times in a given day than those closer to the equator. Table 1 shows the number of images acquired per lake over the study period. Each satellite image was processed using the Color Producing Agents Algorithm (CPA-A) [47] to generate estimates of chlorophyll-a (CHL) concentration, suspended mineral concentration, and colored dissolved organic matter (CDOM) absorption at 443 nm. From these three retrievals, bulk absorption (a) and backscatter (bb) coefficients at 488 nm were computed. The CPA-A uses remote sensing reflectance at the following bands: 443 nm, 488 nm, 531 nm, 547 nm, and 667 nm. Low quality pixels were removed using the following NASA OBPG-provided data quality flags: LAND (pixel over land), HILT (observed radiance very high or saturated), HISATZEN (sensor view zenith angle exceeds threshold), STRAYLIGHT (probable stray light contamination), and CLDICE (probable cloud The Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on board the Aqua satellite was the source of remotely sensed imagery. This imagery was acquired through the National Aeronautics and Space Administration (NASA) Ocean Biology Processing Group (OBPG) (https: //oceancolor.gsfc.nasa.gov/) at Level 2, having been processed by the OBPG using their standard atmospheric correction routines [46]. For each lake, all images collected between 1 January 2003 and 31 December 2018 that intersected any portion of the target were identified and downloaded. Because MODIS Aqua is a polar-orbiting sensor, lakes closer to the poles are observed more times in a given day than those closer to the equator. Table 1 shows the number of images acquired per lake over the study period. Each satellite image was processed using the Color Producing Agents Algorithm (CPA-A) [47] to generate estimates of chlorophyll-a (CHL) concentration, suspended mineral concentration, and colored dissolved organic matter (CDOM) absorption at 443 nm. From these three retrievals, bulk absorption (a) and backscatter (b b ) coefficients at 488 nm were computed. The CPA-A uses remote sensing reflectance at the following bands: 443 nm, 488 nm, 531 nm, 547 nm, and 667 nm. Low quality pixels were removed using the following NASA OBPG-provided data quality flags: LAND (pixel over land), HILT (observed radiance very high or saturated), HISATZEN (sensor view zenith angle exceeds threshold), STRAYLIGHT (probable stray light contamination), and CLDICE (probable cloud or ice contamination). The NASA OBPG CLDICE flag is triggered when infrared reflectance (IR) exceeds a pre-defined threshold which is set to be much greater than that expected from liquid water. While small fragments (relative to the 1 km 2 pixel area) of ice or clouds may contribute to a given pixel's reflectance value, if the IR threshold is not exceeded the pixel value is considered a quality water pixel. Derived products were projected to the World Geodetic Coordinate System of 1984 (WGS84) and saved as GeoTIFF files.
Using methods outlined in Fahnenstiel et al. [38], monthly mean images were created for the CHL, a, and b b parameters. Because the sensor can observe a given lake multiple times within one day due to overlapping swaths, daily mean images were created by averaging overlapping pixels within all images from a given day. Monthly mean images were generated by averaging all daily mean images within the given month. For the CHL product, two additional filters were employed to remove poor retrievals due primarily to invalid reflectance spectra that infrequently passed the quality flag filtering. First, CHL concentrations less than 0.15 mg/m 3 were excluded. Second, daily mean CHL concentrations of more than 7 standard deviations from the monthly mean for a given pixel were excluded, and the monthly mean was re-computed [38]. The proportion of pixels excluded using these filters varied by lake and are shown in Table 2. The remote sensing-derived monthly mean data were used to estimate the light attenuation coefficient in the visible domain (350-700 nm), K PAR , and photic zone depth (the depth where 1% of surface light remains), z 1% , using a methodology outlined in Lee et al. [48]. This approach models light attenuation as a function of the absorption and backscattering coefficients of the water (a, and b b ) and solar zenith angle. Hourly values of K PAR and z 1% were derived using the monthly averaged coefficients and hourly solar zenith angle was determined based on the Julian date and latitude [49]. Latitude was set to be the approximate geometric center of each lake.
Monthly mean lake surface temperature (LST) data were also generated for each lake. For the Laurentian Great Lakes, daily composites from the Great Lakes Surface Environmental Analysis (GLSEA) [50] Version 2 dataset were used to generate monthly averaged LST images. For the other six lakes, monthly mean images were generated using the MODIS Aqua 4 µm Sea Surface Temperature product. All images intersecting each lake were downloaded from the NASA OBPG, having already been processed to Level 2. Several filters were applied to remove low quality data. Pixels failing the NASA-provided SST_CLOUD (pixel failed the cloud alternating decision tree) and SSTRANGE (retrieval outside the physically realistic range) flags were identified and discarded. Monthly average LST data were computed using a weighted average of the data within that month. Weights were assigned based on the NASA-provided quality flag (qual_sst4), with best quality data (level 0) weighted highest and lowest quality data (level 3) weighted lowest. Failed data (level 4) were discarded. Any remaining gaps in the monthly mean image were filled with the lake-wide median temperature and a 9 × 9 Gaussian smoothing filter was applied.
Hourly photosynthetically active radiation (PAR) data for each lake were generated from data acquired through the National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Protection (NCEP). Data from January 2003 through March 2011 were part of the Climate Forecast System Reanalysis (CFSR) dataset [51] and data from April 2011 through December 2018 were part of the Climate Forecast System Version 2 (CFSv2) dataset [52]. CFSR data has a resolution of approximately 31 × 31 km while CFSv2 data has a resolution of approximately 20 × 20 km. The hourly downward shortwave radiation flux product was downloaded for all days in our analysis period for each lake. Shortwave radiation flux was converted to PAR flux with a conversion factor of 0.368 [53] and then to photon flux with a conversion factor of 4.56 [38].
Bathymetric maps for each lake were also acquired or generated from the available data. Laurentian Great Lakes bathymetric maps were acquired from the National Oceanic and Atmospheric Administration (NOAA) [54][55][56][57][58]. The bathymetric map for Lake Victoria was created by Hamilton et al. [59]. Digital bathymetric contour maps were acquired for Lakes Malawi [60], Tanganyika (courtesy of tcarta.com), and Baikal [61]. Contour maps were digitized from published images of Great Slave Lake [62] and Great Bear Lake [63]. Gridded bathymetric maps were created from the contour maps using the "Topo to Raster" tool in ArcMap.
Phytoplankton production was estimated using the Great Lakes Production Model (GLPM) which is a mechanistic approach with technical details described by Lang and Fahnenstiel [64]. This model approach was applied in the Great Lakes using remote sensing data by Fahnenstiel et al., who reported good agreement between GLPM production estimates and field data for Lakes Michigan, Huron and Superior (type II regression model y = 0.98X − 86, R 2 = 0.76, p < 0.001, n = 25) [38]. As implemented in [38], this model uses a combination of parameters that are either remotely sensed or derived from empirical relationships. The remotely sensed parameters include monthly-averaged CHL and LST, and hourly incident (derived from CFSR/CFSv2) and underwater (derived from K PAR ) irradiance as described above. The empirically-derived production parameters include P max and alpha (initial linear slope at low irradiances). The empirical relationship between P max and LST reported for the Laurentian Great Lakes by Fahnenstiel et al. [38] was modified to reflect relationships across the large freshwater lakes of the world. Using the GLPM, production was calculated through the photic zone on an hourly scale assuming vertically uniform chlorophyll-a concentrations and K PAR values. In cases where the derived photic zone depth exceeded the water depth from the bathymetric data, production was calculated to the depth from the bathymetric data layer.
The hourly production data within a given day was summed to estimate daily areal primary production (mg C/m 2 /day). Daily production data within a given month were averaged to generate monthly mean production maps. Finally, the monthly mean production maps were averaged together in order to generate annual mean daily areal primary production maps for each lake. Gaps in monthly production data indicated that remotely sensed data (either CHL or LST) was not available for those pixels in that month due to ice or cloud cover. Due to extended periods of ice and cloud cover in four of the lakes, a subset of months were used when calculating annual mean production maps for Lake Erie (March-December), Lake Baikal (May-November) and Great Bear Lake and Great Slave Lake (June-November). Mean annual maps were also generated for the input parameters, including CHL, K PAR , LST, and incident irradiance (calculated as the mean total daily irradiance over the year). Mean annual carbon fixation was estimated for each lake by multiplying the mean annual lake-wide production estimate by the lake area, then converted to units of TgC/year. Global meteorological data were downloaded from the MERRA-2 (Modern-Era Retrospective Analysis for Research and Applications, Version 2) dataset [65]. Specifically, the surface wind speed product from the Surface Land Forcings dataset was acquired. Monthly and yearly means were generated from the hourly dataset for each lake.
All analyses were performed in the R software (version 3.2.3). Differences in production between lakes were assessed using analysis of variance (ANOVA), with the Tukey HSD test being used to assess pairwise differences. Significance of time-series trends was assessed using ordinary least squares (OLS) regression [66]. A significance level of 0.05 was used for all analyses.

Results
Two of the GLPM input parameters, P max and alpha, were derived empirically rather than from remote sensing directly. A significant relationship was found between P max and LST in the Laurentian Great Lakes [38], allowing for the estimation of P max from remote sensing-derived LST. Because this analysis includes data from freshwater lakes around the globe, a new relationship between P max and LST was derived using observations from large freshwater systems in Africa and North America (n = 319, sources listed in Table 3). No significant difference was found between the two regions (p-value > 0.05) allowing for a universal relationship between P max and LST to be used for all freshwater lakes in this study ( Figure 2). Table 3. Sources for historic alpha and P max data. The mean alpha value from each study is also included, rounded to the nearest integer.   The alpha parameter was also determined by including data from the other large lakes of the world (sources listed in Table 3). Fahnenstiel et al. [38] estimated the alpha value by month based on an observed relationship in the Laurentian Great Lakes. The larger alpha dataset exhibited little variability and no significant trends were observed relative to ancillary data (i.e., month, temperature, lake, latitude, etc.). Therefore, the mean of the alpha dataset, 3.8 mgC/mg Chl/mol photon/m 2 , was used across all calculations.

Lake
Mean annual lake-wide production and carbon fixation were calculated for all 11 lakes for the years 2003-2018. These estimates, along with the minimum and maximum production are displayed Figure 2. Maximum photosynthetic rate (P max ) plotted against lake surface temperature from large freshwater lakes in Africa and North America (bin averaged 1 • C). The exponential relationship (y = 1.012 × e 0.0515x ) was highly significant (p-value < 0.001, R 2 = 0.86).
The alpha parameter was also determined by including data from the other large lakes of the world (sources listed in Table 3). Fahnenstiel et al. [38] estimated the alpha value by month based on an observed relationship in the Laurentian Great Lakes. The larger alpha dataset exhibited little variability and no significant trends were observed relative to ancillary data (i.e., month, temperature, lake, latitude, etc.). Therefore, the mean of the alpha dataset, 3.8 mgC/mg Chl/mol photon/m 2 , was used across all calculations.
Mean annual lake-wide production and carbon fixation were calculated for all 11 lakes for the years 2003-2018. These estimates, along with the minimum and maximum production are displayed for each lake in Table 4, sorted by mean annual production. Yearly lake-wide production estimates for each lake are available in Table S1. Table 4. Mean, minimum, and maximum annual production by lake, sorted by mean annual production (low to high). The ordinary least squares (OLS) slope of the 16-year time series is also included as well as the mean annual carbon fixation estimate. Annual production differed between lakes (p < 0.001), and a post-hoc analysis revealed eight lake pairs with statistically similar production. Great Slave Lake, Lake Huron, Lake Michigan, and Lake Superior all showed a similar level of production (all p > 0.05). Additionally, Great Slave Lake was similar to Lake Baikal (p > 0.05) and Lake Ontario was similar to Lake Malawi (p > 0.05). All other lake pairs were statistically dissimilar.
The 16-year time series allows for the assessment of trends in lake-wide productivity. Of the 11 large lakes, only three exhibited changes in productivity over the entire 16 year time period; productivity was seen to be significantly increasing in Great Bear Lake (slope = 2.664, R 2 = 0.349, p = 0.016) and Great Slave Lake (slope = 4.845, R 2 = 0.3, p = 0.028) and decreasing in Lake Tanganyika (slope = −5.718, R 2 = 0.369, p = 0.013) (Figure 3). Annual trends were also assessed for the model input parameters (CHL, KPAR, LST, and PAR) as well as wind speed (Table 5). Significant CHL declines were observed in Lakes Tanganyika and Victoria. Significant declines in KPAR were observed in Lakes Tanganyika, Victoria, and Michigan. A significant increase in LST was observed in Lake Baikal. Significant increases in PAR were observed in Lakes Baikal and Victoria. Significant declines in wind speed were observed in Lakes Tanganyika and Victoria. Annual trends were also assessed for the model input parameters (CHL, K PAR , LST, and PAR) as well as wind speed (Table 5). Significant CHL declines were observed in Lakes Tanganyika and Victoria. Significant declines in K PAR were observed in Lakes Tanganyika, Victoria, and Michigan. A significant increase in LST was observed in Lake Baikal. Significant increases in PAR were observed in Lakes Baikal and Victoria. Significant declines in wind speed were observed in Lakes Tanganyika and Victoria. Table 5. OLS trend direction (increasing or decreasing) and p-value for each lake's model input parameters and MERRA-2 derived wind speed.

Discussion
This study is the first to compare primary production in the world's largest lakes using a consistent methodology. Although primary production measurements have been collected and reported on in each of the large lakes over the last fifty-plus years ( [23,[29][30][31][32][33][34][35][36][37][38]), differences in measurement techniques among investigators has made it difficult to compare production values between lakes, or even between studies in the same lake [35,36,39]. Thus, while the results from these past studies are useful, they lack the discrimination necessary for a broad comparison among the world's largest lakes. By using consistent methodology across all lakes, this study allows for the first ever comparison of primary production across the world's largest lakes without the need to consider methodological or spatial extent differences. The use of remote sensing also greatly increases the number of observations both in space and time used to determine our estimates. For example, over 15 million CHL and K PAR observations were used for estimating primary production in each of the lakes, with between 16 and 33 million observations in each of the African Great Lakes; 22-33 million in the near-arctic lakes and 22-114 million in the Laurentian Great Lakes.
A wide range of mean annual primary production values was noted across the large lakes, from under 200 mgC/m 2 /day in Great Bear Lake to over 1100 mgC/m 2 /day in Lake Erie. Based on annual production values, the lakes fall into three groups. Great Slave Lake and Lakes Huron, Michigan, and Superior are all statistically similar, as are Great Slave Lake and Lake Baikal. Great Bear Lake is statistically less productive than the rest. These six least-productive lakes, also the northern-most lakes studied, will be considered our oligotrophic lakes. At a higher level of production, Lakes Ontario, Malawi and Tanganyika make up our mesotrophic lake group. While, Lake Tanganyika is statistically more productive than Lakes Ontario and Malawi, all three lakes have mean annual production within 100 mgC/m 2 /day of each other. Finally, Lakes Erie and Victoria are grouped together. Despite being statistically dissimilar, these two lakes have much higher production than the other lakes and each has a history of eutrophication due to human development and land use change [25,27,68]. These two lakes will be considered our eutrophic lakes. Being able to make these between-lake productivity comparisons is an interesting result that is only made possible by using a consistent approach across all lakes.
The use of the GLPM approach for estimating Laurentian Great Lakes production using remote sensing-derived input data was validated by Fahnenstiel et al. [38], who found good agreement between estimates derived from remote sensing and field measurements from three of the Lakes (Michigan, Huron, and Superior). This study made several slight modifications to the methods outlined in [38], including updates to the empirical relationships to reflect the broader study region and a different approach for the calculation of K PAR . Despite these changes, our production estimates agree closely with those reported in [38]. Fahnenstiel et al. reported mean annual production for Lakes Huron, Michigan, and Superior from 2010-2013 as 216, 259, and 228 mgC/m 2 /day, respectively [38]. Over the same time span, this study reported production at 231 mgC/m 2 /day for Lake Huron (7% difference), 270 mgC/m 2 /day for Lake Michigan (4% difference), and 246 mgC/m 2 /day for Lake Superior (8% difference).
Given considerations for methodological differences, overall our production estimates compare favorably to other published estimates in each of the lakes studied. Sterner reported lake-wide annual production for Lake Superior at 257 mgC/m 2 /day for the 2006-2008 time period [37], while our study estimated production at 249 mgC/m 2 /day over the same range. There are few recent estimates of lake-wide annual production in either Lake Ontario or Lake Erie. Lake Ontario production was estimated between 728 and 761 mgC/m 2 /day for May-October 1987-1991 [31]. Averaging across the same six months in our data record, our mean daily production was within that range, at 735 mgC/m 2 /day. More recent estimates of production in Lake Ontario [69,70] were excluded due to significant methodological differences. Recent production estimates for Lake Erie are all reported at the basin scale (i.e., western, central, eastern) and so are not directly comparable with the lake-wide estimates from this study due to the lake's significant between-basin variability [71,72].
Similarly, our estimates for the African Great Lakes and Lake Baikal compare favorably with historical estimates. In Lake Tanganyika, Bergamino et al. [45] used remote sensing to derive an estimate of annual lake-wide production of 646 mgC/m 2 /day in 2003 and Stenuite et al. [36] estimated annual production at 340 mgC/m 2 /day at a northern station and 560 mgC/m 2 /day at a southern station. Our estimate of annual lake-wide production at 591 mgC/m 2 /day is within range of these estimates. Guildford et al. [35] reported annual production in Lake Malawi at 462 mgC/m 2 /day, consistent with this study's estimate of 500 mgC/m 2 /day. Additionally, production in Lake Victoria from 2000-2001 was estimated to be 825 mgC/m 2 /day [33], similar to our estimate of 854 mgC/m 2 /day. Even though there are few recent estimates for Lake Baikal, and these were mostly regional (focused on the southern basin), the historical estimates are within the range of our estimates. Yoshida et al. [32] summarized much of the existing published production data for southern Lake Baikal. Reported values of mean annual production in the region included 205 mgC/m 2 /day in 1999-2000 [32] and 448 mgC/m 2 /day in 1981-1985 [73], bracketing our lake-wide estimate of 327 mgC/m 2 /day.
The one region where our estimates do not agree as well with historical estimates is the near-arctic lakes. Published estimates for Great Bear Lake (11 mgC/m 2 /day [29]) and Great Slave Lake (27 mgC/m 2 /day [29]; 41-110 mgC/m 2 /day [30]) were well below those estimated by this study. Because these estimates were from sampling completed in the 1970s and 1980s with very different methods, it is not surprising our estimates are higher [38].
Despite the significant environmental changes encountered in the world largest lakes in the last twenty years [15,16,74], only three lakes showed significant changes in production over the study period. Only Great Bear and Great Slave lakes and Lake Tanganyika exhibited significant changes in the 2003-2018 period. It appears that these production changes are at least partially due to climate change. This conclusion is supported by both an analysis of the model input parameters and recent documentation of environmental changes. Neither of the near-arctic lakes showed significant changes in any of the input parameters, but LST and PAR (possible climate indicators) increased during the study period in both lakes (p-values ranging from 0.08-0.21). While not significant individually, these combined climate effects appear to have driven the observed production increase. A climate-driven increase in primary production in these near-arctic lakes is not surprising and is consistent with projected increases due to changes in precipitation, temperature and duration of ice cover [17]. NOAA's 2018 Arctic Report Card indicates that these two lakes froze up 2-3 weeks later than average in 2017-2018 and that ice cover duration has been below the historic average in 9 of the prior 14 years [18], further supporting the notion that climate change is impacting these lakes. Even though the recent changes in Lake Tanganyika were opposite (decrease) to those in Great Bear and Slave Lakes, they also appear to be the result of recent climate change. The decline we noted appears to be part of a larger trend dating back to at least to the mid-20th century [12,13,36] which has been linked to a changing climate, with a combination of increasing temperatures and decreasing wind speeds leading to a reduction in mixing events which provide necessary nutrients for phytoplankton production [12,13,75]. The obvious climate change-related input parameters, LST and PAR, did not change significantly in Lake Tanganyika over the 16 year study period, but significant declines were observed in wind speed, CHL and K PAR , consistent with previously reported environmental changes and drivers of production decline.
With their close proximity, it may be surprising that the climate trends impacting Lake Tanganyika did not also lead to a production decline in nearby Lakes Malawi and Victoria. In Lake Malawi, this was likely due to the lack of significant change in any of our climate related parameters (LST, PAR, or wind speed). Lake Malawi also likely responds differently to environmental forcing due to increased seasonality in air temperatures and precipitation [76]. At the annual scale of this study, we were unable to evaluate whether seasonal production was changing, but our results do support the uniqueness of each African large lake. Primary production in Lake Victoria also did not change during the study period, and similar to in Lake Malawi, this was likely related to a different response to environmental variables. As in Lake Tanganyika, CHL, K PAR , and wind speed declined over the study period, though Lake Victoria experienced a significant increase in PAR as well, which could have offset the other changes. Lake Victoria is also closer to the equator and much shallower than Lake Tanganyika, resulting in different mixing dynamics such that we would not expect the two lakes to respond the same way to climate change [74]. Due to a combination of anthropogenic factors, climate change, and invasive species, Lake Victoria has a long history of eutrophication extending to the early 20th Century [27,68,77], but it was a period of low wind stress and increasing temperatures beginning in the 1970s that transformed the ecosystem to being more productive [28]. Because production did not change and chlorophyll-a concentrations declined over the last 18 years the eutrophication period in Lake Victoria may be over.
The lack of any significant change in productivity and almost all input parameters for the Laurentian Great Lakes (the only parameter that showed significant change was K PAR in Lake Michigan) is surprising as many previous investigators have documented significant changes in production and climate related parameters over the past two decades. For example, significant summer LST warming trends have been identified for many of the Laurentian Great Lakes [14,16,78] and a decline in ice cover duration has been reported in all of the lakes [14,16]. It should be noted that LST may not be the best indicator of climate warming since near surface temperatures can at times have opposing trends to total lake heat content, which is a better indicator of climate warming [79]. Moreover, significant spatial variability has been noted in LST and ice cover trends within lakes [16], and thus, our annual lake-wide trends in production and input parameters may miss some changes on finer spatial and temporal scales. However, because of the accuracy and precision of many of our input parameters, i.e., LST, the lack of significant trend in nearly all input parameters, including climate indicators (LST, PAR, wind speed), is noteworthy.
Invasive species have also been noted to impact production and phytoplankton biomass in the Laurentian Great Lakes. Specifically, the introduction of dreissenid mussels (zebra mussels, Dreissena polymorpha, and quagga mussels, Dreissena bugensis) has led to many significant changes throughout the lakes. First identified in the lower Laurentian Great Lakes in the late-1980s [20], these mussels quickly established themselves throughout the lakes to become the dominant benthic organisms by the late-1990s to mid-2000s [22,80,81]. Lake Superior is the only Laurentian Great Lake that has not suffered a widespread invasion of these mussels [82,83]. Significant declines in phytoplankton biomass and production through the 1990s were linked to dreissenid filtering both in Lake Erie [22,84,85] and Lake Ontario [21,70]. With effects extending into our study period, the mussels have also been deemed responsible for reducing spring and winter phytoplankton blooms in Lake Michigan between the 1990s and late-2000s [23,86] and increasing water clarity in Lakes Huron and Michigan from 1998 through 2010 [24]. Water clarity in Lake Michigan was observed to be significantly increasing across the study period, but no other significant changes were noted in CHL or K PAR for the other lakes. While the impact of these mussels is clear from the literature, our study period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) primarily captures the post-invasion phase when the quagga mussel had already been established as the dominant benthic organism and many of the effects were already manifested. While the mussels are likely still impacting the ecosystem, it is possible that the effects of their continued filtration are now serving to offset production changes related to other anthropogenic factors (i.e., climate change, eutrophication, etc.).
Finally, the lack of our ability to resolve trends in all of the lakes does not mean that actual trends are not present, but that the approach used in this study may limit our ability to resolve small changes. For example, there was significant interannual variability in the annual production estimates for all lakes, with coefficient of variation ranging from 5-14%. Lake Baikal had the largest coefficient of variation of all lakes that did not exhibit a significant production trend (11.2%), indicating that it would be difficult to identify production changes of less than 20%. Additionally, the relatively limited time span of the remote sensing data used may be too short to capture environmental trends that manifest over decades. For example, an analysis of changes in δ 13 C org levels in Lake Superior sediments revealed that production has been increasing over the last century [87]. Extending the study period by adding subsequent years may help to resolve these issues and determine anthropogenic influences. Another factor that may contribute to our lack of trends in production is the potentially offsetting interaction of input parameters (Table 5). For example, no significant trends were observed in the mean annual production for Lakes Baikal and Ontario. However, when looking at the annually-averaged input parameters, both lakes showed significant or near-significant increases in climate-impacted parameters LST and PAR with offsetting decreases in CHL and K PAR . Thus, climate change may be impacting Lakes Baikal and Ontario, but we were not able to determine any significant change in production due to possible impacts of invasive species and/or eutrophication. On the other hand, Lake Huron also showed no significant production trend, and an inspection of the input parameters also revealed no significant changes.
This study is the first to compare annual lake-wide primary production in the world's largest lakes using the same methodology. Some of these lakes have been regularly studied, while others had not been measured in over 30 years. This comparison would not have been possible without remote sensing, which eliminates much of the temporal and financial barriers to widespread sampling. However, field measurements will continue to be needed in order to validate the remotely sensed retrievals and to improve the algorithm. For instance, the P max and alpha measurements used were between eight and twenty-three years old and only represented six of the eleven lakes studied. Collecting a more representative set of these measurements (across all lakes and seasons) could lead to an improved set of estimates as better relationships are developed for photosynthetic parameters and environmental variables. Additionally, measurements of inherent optical properties in the lakes could help to improve the accuracy of the CHL and K PAR input parameters [88]. Beyond this study, further exploration of the dataset will serve to better characterize the lakes and their productivity trends. Many studies have noted that production in large lakes is spatially and seasonally dependent [23,35,45], and these previously reported trends could be re-evaluated with our more expansive data set or new relationships could be identified. Focusing on specific regions or seasons could also help to reduce the uncertainty in our estimates, allowing for easier evaluation and explanation of trends.

Conclusions
This study used a satellite remote sensing based approach to estimate phytoplankton primary production in 11 of the world's largest lakes from 2003-2018. Trends in annual production were analyzed and revealed only three of the eleven lakes experienced statistically significant changes over the 16 year period with Great Bear and Great Slave Lakes exhibiting significant increases in production while Lake Tanganyika experienced a significant decline. These documented trends were linked to changes in climate, most notably increases in water temperatures and solar radiation and a reduction in wind speed. The absence of climate change impacts on the other eight lakes may be due to the limited temporal observation period (16 years) as well as the offsetting interactions of other production model input variables (e.g., decreasing CHL biomass and increasing water clarity). While we focused on climate induced changes in production, further investigation is needed to explore linkages between changes in other local and regional environmental factors (e.g., land-use/land-cover) and production patterns in large lakes. Finally, remote sensing time-series assessments of primary production for additional freshwater systems representing a larger range in size and trophic status would be extremely insightful toward exploring the impacts of a changing climate on a wide range freshwater ecosystems.