Spatio-Temporal Variability of Phytoplankton Primary Production in Baltic Lakes Using Sentinel-3 OLCI Data

: Phytoplankton primary production (PP) in lakes play an important role in the global carbon cycle. However, monitoring the PP in lakes with traditional complicated and costly in situ sampling methods are impossible due to the large number of lakes worldwide (estimated to be 117 million lakes). In this study, bio-optical modelling and remote sensing data (Sentinel-3 Ocean and Land Colour Instrument) was combined to investigate the spatial and temporal variation of PP in four Baltic lakes during 2018. The model used has three input parameters: concentration of chlorophyll- a , the di ﬀ use attenuation coe ﬃ cient, and incident downwelling irradiance. The largest of our studied lakes, V õ rtsjärv (270 km 2 ), had the highest total yearly estimated production (61 Gg C y − 1 ) compared to the smaller lakes Lubans (18 Gg C y − 1 ) and Razna (7 Gg C y − 1 ). However, the most productive was the smallest studied, Lake Burtnieks (40.2 km 2 ); although the total yearly production was 13 Gg C y − 1 , the daily average areal production was 910 mg C m − 2 d − 1 in 2018. Even if lake size plays a signiﬁcant role in the total PP of the lake, the abundance of small and medium-sized lakes would sum up to a signiﬁcant contribution of carbon ﬁxation. Our method is applicable to larger regions to monitor the spatial and temporal variability of lake PP.


Introduction
Lakes and reservoirs play a key role in global carbon cycling, either as carbon sinks [1] or sources [2][3][4]. According to Verpoorter et al. [5], there are about 117 million lakes larger than 0.002 km 2 , covering 3.7% of the Earth's non-glaciated land surface. Despite their relatively small surface area, compared to the open ocean, lakes make a significant contribution to the global carbon budget, being extremely active sites for terrestrial carbon transport, transformation, and storage [6]. Understanding the drivers of carbon cycling and the ways climate change affects them are important goals of ecology and biogeochemistry in the present day [7].
Phytoplankton primary production (PP) is the process through which algae and bacteria fix inorganic carbon and transform it into organic matter; it is one of the ways in which inorganic carbon enters the aquatic ecosystem. This forms the basis of aquatic food webs, generates food sources for

Study Sites
The study area included three Latvian and one Estonian lake, all of which are located in the northern boreal zone. The locations of the lakes are shown in Figure 1, and their geographical characteristics are shown in Table 1.  Lake Razna is located in the Latgale highlands in Latvia and has a surface area of 57.6 km 2 and catchment area of 229 km 2 . It is the country's second-largest lake by surface area and the largest by water volume (0.46 km 3 ). According to the national monitoring database, Lake Razna is classified as a shallow lake with clear water [53], although the maximum depth reaches 17 m (mean depth of 7.0 m). The lake is mesotrophic; the Secchi depth (SD) was measured from 3.2 to 6.6 m (median 5.5 m) and Chl a from 1.4 to 10.9 mg m −3 (median of 6.4 mg m −3 ) during April to October 2018 [37].
Lake Lubans is situated in the center of the Eastern Latvian Lowland. The catchment area of the lake is 2040 km 2 . At a water level of 90 m above sea level, the surface area of the lake is 80.7 km 2 . However, the surface area of the lake can fluctuate from 25 to 100 km 2 due to an artificially regulated water level [54]. Lake Lubans is classified as a very shallow lake [53], where the maximum depth is 3.5 m and a the mean depth is 1.6 m. The lake is eutrophic and macrophyte-dominated; the SD was measured at 0.4-1.2 m (median of 0.9 m) and Chl a at 8.6-63.1 mg m −3 (median of 16.8 mg m −3 ) during April to October 2018 [37].
The large Lake Võrtsjärv (surface area of 270 km 2 ) is in the central part of Estonia. It is the country s second-largest lake, with a catchment area of 3374 km 2 . Lake Võrtsjärv is also a shallow lake, with a maximum depth of 6.0 m (mean depth of 2.8 m). The unregulated water level can fluctuate annually up to 3.2 m (mean 1.4 m) [55]. In this eutrophic lake, the SD varied between 0.5 and 1.0 m (median of 0.6 m) and Chl a varied between 22.2 and 63.8 mg m −3 (median 33.1 mg m −3 ) during April to October 2018 [37].
Lake Burtnieks is situated in the northern part of Latvia. It is the fourth-largest lake in Latvia, with a surface area of 40.2 km 2 and a catchment area of 2215 km 2 . Lake Burtnieks is a shallow eutrophic lake with a maximum depth of 4.3 m (mean depth of 2.4 m). The lake has brown water, where the SD was measured 0.4-1.2 m (median of 1 m) and Chl a at 6.7-117.1 mg m −3 (median of 19.1 mg m −3 ) during April to October 2018 [37].

Satellite Data
Cloud-free or partially cloud-free OLCI Level-1 full resolution images from onboard the Sentinel-3A and -3B satellites in 2018 were used in this study. Images were downloaded from the Copernicus Online Data Access platform [56] and processed with the atmospheric correction processor Case 2 Regional Coast Color (C2RCC) (v1) [57] using Sentinel Application Platform (SNAP) (v6) software developed by Brockmann Consult, SkyWatch, and C-S [58].
The retrieval of the water-leaving reflectance from OLCI Level-1 images involved three-step processing: (1) subsetting the case study lake from the entire scene using shape files of each lake; (2) applying the multisensor pixel identification tool (IdePix) and using in further analysis only the cloud-free inland water pixels by excluding pixels flagged as IDEPIX_CLOUD, IDEPIX_CLOUD_AMBIGUOUS, IDEPIX_CLOUD_SURE, IDEPIX_CLOUD_BUFFER, or IDEPIX_CLOUD_SHADOW; and (3) applying a C2RCC atmospheric correction processor using default parameters (except for temperature, which was set as 15 • C, and for salinity, which was set as 0.001 PSU) and excluding from further analysis pixels that were flagged as Kd489_OOR (K d,489 product is out of range) or Cloud_risk. Only the images with ≥50% of the valid lake pixels were included in the further analysis.
Since standard remote sensing Chl a and K d products often fail [59][60][61][62] in these optically complex inland waters, an alternative solution suggested for this region was used (detailed explanation are in Sections 2.4.1 and 2.4.2).

The Primary Production (PP) Model
The PP model, originally developed by Arst et al. [20,52] and further adapted for remote sensing data by Kauer et al. [41] and Soomets et al. [42] was used in this study. The main principle of this model is that the PP is a function of photosynthetically absorbed radiation and the quantum yield of carbon fixation [63]: where PP(z) is production in mg C m −3 h −1 at depth z, Ψ is the factor 12 000 for converting moles of carbon to milligrams of carbon, Q * PAR (z) is photosynthetically absorbed radiation at depth z determined on the basis of scalar quantum irradiance (in mol photons m −3 h −1 ), and F PAR (z) is the quantum yield of carbon fixation (mol C (mol photons) −1 ) in the PAR at depth z. The PP model has three input parameters: (1) Chl a in mg m −3 , (2) K d,PAR in m −1 , and (3) q PAR in mol photons m −2 h −1 . In this study, we followed the calculation scheme of the PP model modified for remote sensing data by Soomets et al. [42], and we used the integrated values of primary production (PP int , mg C m −2 h −1 ) by integrating PP(z) over the mean depth of the lake for all valid pixels.
Based on the PP int of the model, the following further calculations were made: • The total lake PP int (PP lake , mg C h −1 ) for the entire lake. As all available images with ≥50% of the valid lake pixels were used in the study, PP lake was calculated by summarizing the available pixel values of PP int , and then the average of the available pixels of a given day was taken to substitute all of the missing surface-area pixels.

•
The average areal PP int of the lake (PP aver , mg C m −2 h −1 ) was taken as the average of the PP lake over the surface area of the lake. For this parameter we also calculated the standard deviation (±) over the study period (April-October).

•
The daily PP lake (PP lake,day, mg C d −1 ) was calculated by multiplying PP lake , the photoperiod (in h) of a given day, and the coefficient of 0.75 to take into account a daily light curve.

•
The average daily areal PP int for the entire lake (PP day , mg C m −2 d −1 ) was taken as the average of the PP lake,day over the surface area of the lake.

•
The monthly average PP day (PP day,aver , mg C m −2 d −1 ) was taken as the average of the PP day during one month.

•
The monthly productivity estimates of the entire lake (PP month , Gg C month −1 ) were taken as the sums of PP lake,day for each month during the study period. If the PP lake,day value was missing due to the absence of overpass of the satellite or cloudy conditions, then inter-or extrapolation was used to calculate missing the PP lake,day value.

•
The annual productivity estimates of the entire lake (PP year , Gg C year −1 ) were taken as the sums of PP month during April-October 2018. We assumed that the PP during the rest of the months during 2018 is zero or close to zero due to ice and snow cover and not significant for the annual estimation of PP. Besides, there were no cloud-, ice-, or snow-free images available for the studied area from January to March, and November to December 2018.

•
The annual average daily areal productivity (PP year,aver , mg C m −2 d −1 ) was calculated by dividing PP year by the number of days in a year and the surface area of the lake.

Input Parameters of the PP Model
The PP model has three input parameters: Chl a, K d,PAR , and q PAR . To study the seasonal and spatial variatios of PP in lakes, we used OLCI data for Chl a and K d,PAR and in situ data for q PAR for the PP model. The details of each input are described below.

Chlorophyll-a Concentration (Chl a)
For Chl a, the OWT guided approach by Soomets et al. [37] was used to retrieve Chl a values from OLCI water-leaving reflectance spectra. Classification of OWTs by Uudeberg et al. [64] divides boreal region inland and coastal water into five OWTs: (1) Clear OWT corresponds to waters with low OSC concentrations and high water transparency; (2) Moderate OWT corresponds to slightly higher OSC concentrations, but none of them dominate; (3) Turbid OWT corresponds to waters where TSM dominates; (4) Very Turbid OWT corresponds to waters where Chl a dominates, and it is associated with phytoplankton blooms; and (5) Brown OWT corresponds to waters where the CDOM dominates. The pixel was marked as "Unclassified" and excluded from further analyses if the reflectance spectrum of the pixel was "abnormal" in the blue part due to strong influence from sources other than water.
The OWT guided approach to retrieve Chl a values means that firstly, OWT is assigned to every valid OLCI image water pixel using a method by Uudeberg et al. [64]; and secondly, the most suitable algorithm developed by Soomets et al. [65] for each OWT is applied according to the pixel's OWT. The Chl a algorithm used for each OWT is shown in Equation (2): where R TOA,x is the OLCI top-of-atmosphere reflectance at the specific band that is labeled with its central wavelength.

The Underwater Light Diffuse Attenuation Coefficient (K d,PAR )
For K d,PAR , the calculations from OLCI image data consisted of two main steps. Firstly, the underwater light diffuse attenuation coefficient at 490 nm (K d,490 ) was calculated from water-leaving reflectance spectra using weighted band ratio algorithms by Alikas et al. [66]: where W is a calculated weight that is used to combine the two separate K d,490 algorithms. It is calculated according to Equation (4): where R 560 and R 709 are the OLCI water-leaving reflectance values on specific bands. The two K d,490 algorithms from Equation (3) are calculated using Equation (5): where R λ1 and R λ2 mark the water-leaving reflectance values on specific OLCI bands (490, 560, or 709 nm). K w,490 is the pure water light attenuation coefficient at 490 nm and is 0.016 m −1 ; A and B are coefficients derived by linear regression analysis for the band ratio algorithms: −0.9461 and 0.4305 for K d,490,( 490 709 ) , and −1.1198 and 1.4141 for K d,490,( 560 709 ) . Secondly, we converted K d,490 to K d,PAR using Equation (6):

Incident Planar Downwelling Irradiance (q PAR )
The hourly average diffused sky irradiance for PAR (in µmol m −2 s −1 ) was measured in situ with the SkySpec spectroradiometer [67] at the Station for Measuring Ecosystem-Atmosphere Relations, SMEAR-Estonia, a research station in Järvselja, Estonia (58 • 16 40.60 N, 27 • 18 30.84 E), during 2018. For the PP model q PAR input, we used the measured daily maximum values converted to mol photons m −2 h −1 for each day we had OLCI data. We used the same maximum daily q PAR values for all of the studied lakes because they are located close enough to each other (the maximum distance from the station to the lake was about 220 km).

Spatial Variability of PP
The spatial variability of PP int during 2018 was studied via derived PP int spatial variability maps from OLCI cloud-free or partially cloud-free images; altogether (for all four lakes), we generated 161 PP int spatial maps (approximately 40 for each studied lake). The used PP model has three input parameters, of which two were spatiality variables: Chl a and K d,PAR . Since Chl a was determined using the OWT guided approach, the examples of the resulting PP int (mg C m −2 h −1 ) together with the derived inputs of the PP model, K d,PAR , Chl a, and OWT, are shown in Figures 2-5.

Lake Razna
The clear water Lake Razna showed little spatial variability in its PP int . In spring and summer, the likelihood of higher PP int was at the northern shoreline of the lake ( Figure 2). However, slightly higher productivity in the southern part was often present in autumn. The Clear OWT was very dominant in Lake Razna, and a majority of the pixels had low Chl a (mean 11.1 ± 17.3 mg m −3 ) and K d,PAR (mean 0.62 ± 0.57 m −1 ) values; therefore, the productivity of the lake was also low. The PP aver over the study period was 51.2 ± 14.9 mg C m −2 h −1 .

Lake Lubans
The eutrophic Lake Lubans showed high spatial variability in its bio-optical characteristics; the waters were classified mostly Turbid and Very Turbid throughout the study period. A very distinguishing pattern for Lake Lubans was the decrease of the open water area during vegetation season due to the growing vegetation. This was clearly seen in late August and September, where the open-water area was diminished markedly (Figure 3, 19 September 2018). The mean Chl a was 42.4 ± 40.0 mg m −3 , with a much lower median value (29.4 mg m −3 ), and the mean K d,PAR was 2.52 ± 1.08 m −1 . The PP aver over the study period was 90.9 ± 53.9 mg C m −2 h −1 . 14.04.2018

Lake Võrtsjärv
The largest of the studied Lake Võrtsjärv also showed a relatively high spatial variability of PP int . The PP int was often higher in the southern part of the lake, and in some cases, the inflow of a river to the lake was clearly visible in the southwestern part, where the PP int values are slightly lower than in the surrounding waters ( Figure 4, 19 September 2018). The dominant OWT was either Turbid or Very Turbid. The mean Chl a was 33.4 ± 20.4 mg m −3 but was much higher in the bloom condition (Figure 4), and the mean K d,PAR was 2.15 ± 0.55 m −1 . The PP aver over the study period was similar to Lake Lubans: 100.0 ± 65.1 mg C m −2 h −1 .

Lake Burtnieks
The smallest and the brownest water lake of the current study showed the highest variability in all of the bio-optical characteristics and PP int . Higher productivity of the entire lake was shown especially in late summer. Often, there is a distinct shoreline of higher productivity around the lake ( Figure 5). The dominant OWT was either Turbid, Very Turbid, or Brown, but the Moderate OWT was also frequently present. The mean Chl a was 43.5 ± 37.3 mg m −3 , with a much lower median value (27.7 mg m −3 ), and the mean K d,PAR was 2.58 ± 1.02 m −1 . The PP aver over the study period was 139.2 ± 108.6 mg C m −2 h −1 . 14.04.2018

Temporal Variability of the PP
We calculated the average daily areal production (PP day ) and the total daily production of the lake (PP lake,day ) from the OLCI data for the entire lake. For the dates that we did not have OLCI data, we interpolated (or extrapolated to April 1 and October 31 to get full months) the modelled PP lake,day values. In the period of April to October, we had about 40 days of satellite data (Razna 38, Burtnieks 40, Lubans 42, Võrtsjärv 41), but the satellite data was not distributed equally. The highest number of cloud-free images was acquired in May, whereas in September and October, we had only a few images to use. From November to April, we either did not have any cloud-free images or the lakes were ice-covered. The temporal variability of the PP day shows different patterns for each lake ( Figure 6). In Lake Võrtsjärv and Lake Burtnieks, the PPday showed a clear two-peak pattern over the vegetation period. The first productivity peak was in late spring/early summer (first half of June), and the second productivity rise was in late summer (end of August and September) ( Figure 6). The main difference in the seasonal pattern of those lakes was the duration of the second production peak. In Lake Võrtsjärv, the late summer peak lasted longer (1.5 months) and was not as intense as the spring one. In Lake Burtnieks, both productivity peaks lasted about 3 weeks and the second peak was more intense (Figure 6).
The temporal variability of the PPday had rather similar patterns in Lake Razna and Lubans: both of those lakes are missing the summer "clear" period between the two high productivity periods. Although in Lake Razna, the maximum productivity occurred in early summer, while in Lake Lubans it was in late summer ( Figure 6).
Another very characteristic feature of Lake Burtnieks was the large range of PPday values. For example, the range of PPday was over 4600 mg C m −2 d −1 for the study period, being 4817 mg C m −2 day −1 on 23/08/2018 and 156 mg C m −2 d −1 on 15/10/2018. Other lakes showed much smaller ranges of PPday, with Lake Razna having the smallest range of 950 mg C m −2 d −1 . Lake Razna also had the lowest maximum PPday values (1177 mg C m −2 d −1 ), while in Lake Burtnieks it was 4 times higher. To compare the monthly productivity of the studied lakes, we calculated the monthly averages for PPday (PPday,aver) (Figure 7a). Here, as before, Lake Burtnieks had the highest production rates, but only for the three peak-productivity months (June and August-September). The temporal variability of the PPday,aver in different lakes was demonstrated again, with Lake Lubans having the highest PPday,aver in July, when Lake Võrtsjärv and Lake Burtnieks had a lower PPday,aver (Figure 7a). In Lake Võrtsjärv and Lake Burtnieks, the PP day showed a clear two-peak pattern over the vegetation period. The first productivity peak was in late spring/early summer (first half of June), and the second productivity rise was in late summer (end of August and September) ( Figure 6). The main difference in the seasonal pattern of those lakes was the duration of the second production peak. In Lake Võrtsjärv, the late summer peak lasted longer (1.5 months) and was not as intense as the spring one. In Lake Burtnieks, both productivity peaks lasted about 3 weeks and the second peak was more intense (Figure 6).
The temporal variability of the PP day had rather similar patterns in Lake Razna and Lubans: both of those lakes are missing the summer "clear" period between the two high productivity periods. Although in Lake Razna, the maximum productivity occurred in early summer, while in Lake Lubans it was in late summer ( Figure 6).
Another very characteristic feature of Lake Burtnieks was the large range of PP day values. For example, the range of PP day was over 4600 mg C m −2 d −1 for the study period, being 4817 mg C m −2 day −1 on 23/08/2018 and 156 mg C m −2 d −1 on 15/10/2018. Other lakes showed much smaller ranges of PP day , with Lake Razna having the smallest range of 950 mg C m −2 d −1 . Lake Razna also had the lowest maximum PP day values (1177 mg C m −2 d −1 ), while in Lake Burtnieks it was 4 times higher.
To compare the monthly productivity of the studied lakes, we calculated the monthly averages for PP day (PP day,aver ) ( Figure 7a). Here, as before, Lake Burtnieks had the highest production rates, but only for the three peak-productivity months (June and August-September). The temporal variability of the PP day,aver in different lakes was demonstrated again, with Lake Lubans having the highest PP day,aver in July, when Lake Võrtsjärv and Lake Burtnieks had a lower PP day,aver (Figure 7a). Another estimate that we compared between the lakes was the total monthly productivity for the entire lake (PPmonth). The PPmonth showed vast differences in lakes. It varied from 0.46 Gg C month −1 (Burtnieks, October) to 14.4 Gg C month −1 (Võrtsjärv, September) ( Figure 7b). Here, the impact of lake size on the overall productivity is clearly shown. For example, in August, even though Lake Burtnieks showed twice the production rate per area unit compared with Lake Võrtsjärv (3137 and 1411 mg C m −2 d −1 , respectively), the total production of Lake Burtnieks for this month was 3.9 Gg C month −1 whilst for Lake Võrtsjärv, which is more than six times larger, it was three times higher (11.8 Gg C month −1 ) (Figure 7). Lake size plays an important role also in the yearly lake productivity. For the yearly production, we added the monthly sums of lake production from April to October. We assumed that the production from November to March was not significant [68]. The yearly production of 2018 for the entire lake varied from 7 Gg C y −1 (Razna) to 61 Gg C y −1 (Võrtsjärv) (Figure 8). Another estimate that we compared between the lakes was the total monthly productivity for the entire lake (PP month ). The PP month showed vast differences in lakes. It varied from 0.46 Gg C month −1 (Burtnieks, October) to 14.4 Gg C month −1 (Võrtsjärv, September) ( Figure 7b). Here, the impact of lake size on the overall productivity is clearly shown. For example, in August, even though Lake Burtnieks showed twice the production rate per area unit compared with Lake Võrtsjärv (3137 and 1411 mg C m −2 d −1 , respectively), the total production of Lake Burtnieks for this month was 3.9 Gg C month −1 whilst for Lake Võrtsjärv, which is more than six times larger, it was three times higher (11.8 Gg C month −1 ) (Figure 7). Lake size plays an important role also in the yearly lake productivity. For the yearly production, we added the monthly sums of lake production from April to October. We assumed that the production from November to March was not significant [68]. The yearly production of 2018 for the entire lake varied from 7 Gg C y −1 (Razna) to 61 Gg C y −1 (Võrtsjärv) (Figure 8).
As we assumed the productivity to be zero or close to that for 5 months of the year, the annual average daily areal productivity (PP year,aver ) in Lake Razna would be 333.0 mg C m −2 d −1 , in Lake Burtnieks, 887.3 mg C m −2 d −1 , in Lake Lubans, 610.2 mg C m −2 d −1 , and in Lake Võrtsjärv, 622.4 mg C m −2 d −1 . This shows that the most productive was the smallest of our studied lakes. Lake Võrtsjärv and Lake Lubans show very similar PP year,aver , although the total productivity of the lakes was very different. As we assumed the productivity to be zero or close to that for 5 months of the year, the annual average daily areal productivity (PPyear,aver) in Lake Razna would be 333.0 mg C m −2 d −1 , in Lake Burtnieks, 887.3 mg C m −2 d −1 , in Lake Lubans, 610.2 mg C m −2 d −1 , and in Lake Võrtsjärv, 622.4 mg C m −2 d −1 . This shows that the most productive was the smallest of our studied lakes. Lake Võrtsjärv and Lake Lubans show very similar PPyear,aver, although the total productivity of the lakes was very different.

Discussion
To estimate the spatial and temporal variability of PP, we used 161 cloud-free OLCI scenes over four different Baltic lakes during 2018. The OLCI images were processed with a C2RCC atmospheric processor, and to determine the PP model inputs Chl a and Kd,PAR from the reflectance spectra of each pixel, the OWT guided approach for Chl a and the weighted band ratio algorithms method for Kd,PAR was used. The PPday was calculated to estimate the temporal variability. Finally, daily, monthly, and yearly PPint estimates for each lake were determined.
The spatial variability of PPint in Lake Razna showed the presence of phytoplankton blooms on the northern shores of the lake (Figure 2) during spring; this pattern is also shown as the PPday peaks in the temporal variability ( Figure 6). In the autumn, the PPint had higher values in the southeastern part compared to the rest of the lake on some days, likely because of the winds that would drive the phytoplankton to the shoreline area. The PPday,aver values of this mesotrophic lake are quite similar to the other studied (eutrophic) lakes in April (Figure 7a) because Lake Razna is significantly deeper than the other lakes; therefore, the PPint values, which are calculated by integration over the mean depth of the lake, are also higher.
Although Kauer et al. [69] demonstrated that the PPday values can significantly vary from day to day due to the changing climate condition; the rapid fluctuations of PPday in the current work are mostly caused by the different number of pixels used for a specific day. As we had the requirement that at least 50% of the lake pixels should be valid, it sometimes resulted that the parts of the lake with higher productivity were filtered out (due to the clouds) and the clear water pixels with low PP values were used to estimate PPday. This phenomenon was more pronounced in Lake Razna and Lake Lubans, where flagging out part of the pixels had larger impact on the PPday estimates due to the PPint spatial variation. In addition to the missing pixels, there were errors in the atmospheric correction when thin cirrus clouds or other fine atmospheric particles were not correctly flagged out by used processors.
In this study, the assumption was made that the lake-bottom influence is negligible in our PP estimations, and no additional corrections were applied to remove the bottom influence from the remote sensing signal due to the lakes water mainly being optically deep waters. In Lake Lubans, Burtnieks, and Võrtsjärv, the water transparency (SD ≤ 1.2 m) is lower than the average depth (1.6 m,

Discussion
To estimate the spatial and temporal variability of PP, we used 161 cloud-free OLCI scenes over four different Baltic lakes during 2018. The OLCI images were processed with a C2RCC atmospheric processor, and to determine the PP model inputs Chl a and K d,PAR from the reflectance spectra of each pixel, the OWT guided approach for Chl a and the weighted band ratio algorithms method for K d,PAR was used. The PP day was calculated to estimate the temporal variability. Finally, daily, monthly, and yearly PP int estimates for each lake were determined.
The spatial variability of PP int in Lake Razna showed the presence of phytoplankton blooms on the northern shores of the lake (Figure 2) during spring; this pattern is also shown as the PP day peaks in the temporal variability ( Figure 6). In the autumn, the PP int had higher values in the southeastern part compared to the rest of the lake on some days, likely because of the winds that would drive the phytoplankton to the shoreline area. The PP day,aver values of this mesotrophic lake are quite similar to the other studied (eutrophic) lakes in April (Figure 7a) because Lake Razna is significantly deeper than the other lakes; therefore, the PP int values, which are calculated by integration over the mean depth of the lake, are also higher.
Although Kauer et al. [69] demonstrated that the PP day values can significantly vary from day to day due to the changing climate condition; the rapid fluctuations of PP day in the current work are mostly caused by the different number of pixels used for a specific day. As we had the requirement that at least 50% of the lake pixels should be valid, it sometimes resulted that the parts of the lake with higher productivity were filtered out (due to the clouds) and the clear water pixels with low PP values were used to estimate PP day . This phenomenon was more pronounced in Lake Razna and Lake Lubans, where flagging out part of the pixels had larger impact on the PP day estimates due to the PP int spatial variation. In addition to the missing pixels, there were errors in the atmospheric correction when thin cirrus clouds or other fine atmospheric particles were not correctly flagged out by used processors.
In this study, the assumption was made that the lake-bottom influence is negligible in our PP estimations, and no additional corrections were applied to remove the bottom influence from the remote sensing signal due to the lakes water mainly being optically deep waters. In Lake Lubans, Burtnieks, and Võrtsjärv, the water transparency (SD ≤ 1.2 m) is lower than the average depth (1.6 m, 2.4 m, and 2.8 m, respectively; Table 1); therefore, lake bottoms were not visible from the water surface and were not influencing the reflectance spectra. A lake's bottom can be visible up to few meters from the shore; however, the OLCI pixel size is 300 m. In the clear water Lake Razna, the situation was handled more carefully. Since water transparency in Lake Razna was high (SD was up to 6.6 m) and the sandy bottom was visible in the east side of the lake (depth up to 6 m), we investigated the bottom's influence before making this assumption. The reflectance spectra retrieved in lake areas where the bottom was visible acted similarly to that described by Vahtmäe et al. [70]. Nevertheless, in our case, the bottom influence was not strong enough to change the optical water type or the retrieval of Chl a.
However, the bottom influence was slightly shown in the retrieved K d product. But since the PP model used in this study is more sensitive to Chl a than K d in clear water cases, the bottom influence is not detectable in the PP int results.
Overall, the PP day results were quite high for Lake Lubans throughout the summer, without having clear PP day peaks as in Lake Burtnieks ( Figure 6). Due to extremely high PP day peaks in Lake Burtnieks, the PP aver was also about 30% higher than in Lake Võrtsjärv and Lake Lubans. The higher PP int values in the shore areas are quite usual in all of the lakes, but the most prominent is in Lake Burtnieks, caused by the shallower, more nutrient-rich waters. Another possible reason for higher values in the shore area is an adjacency effect-bias in the spectral signal due to perturbations by the radiance reflected by the land, which are propagated by the atmosphere into the satellite sensor [71].
In the case of Lake Burtnieks and Lake Lubans, Clear OWT is shown along the shorelines of the lakes, which is likely the result of the atmospheric correction issues due to processing with C2RCC [64,72]. For example, the atmospheric correction has been proven to be complicated in the case of dark water (Brown OWT), with the blue part of the spectrum being strongly overestimated. This can cause false classification to the Clear OWT [64,73]. Another cause could be the effect of the in-water vegetation on the reflectance spectra that leads to misclassification of the pixel. The insufficient quality flagging of the C2RCC and/or poor clear water pixel masking by the IdePix tool can also cause a situation where these pixels are falsely classified as clear water pixels [3]. However, the added extra quality flag (Kd489_OOR) seems to exclude most of the pixels affected by the in-water vegetation for Lake Lubans. Thus, the two lakes have a different source for the OWT misclassification: for Lake Lubans, the cause is the in-water vegetation, and for Lake Burtnieks, the brown water. Despite the issues addressed above, our OWT classification results for 2018 generally agree with the previous study for 2017 [72]. Lake Võrtsjärv typically had higher PP int in the southern part of the lake (Figure 4), which is narrower and has different optical properties than the open northern part. Nõges et al. [68] used modeling to reconstruct a long time-series of the areal productivity in Lake Võrtsjärv. They obtained 10% lower PP year,aver results than the current study: 558 mg C m −2 d −1 . One of the probable causes of this difference could be the different lengths of the study periods: the current study was focused only on the year 2018, but the result of Nõges et al. [68] is the mean productivity value over 27 years. Also, our results regarding the PP day,aver (Figure 7a) largely agree with Nõges et al. [68]. Their Figure 4 shows similar monthly values, except for September, where the PP day,aver is almost three times lower. This difference might be explained by the yearly differences in weather conditions. Although the September PP day,aver values were mostly around 600-1000 mg C m −2 d −1 for the years 1982-2009, it was approximately 2000 mg C m −2 d −1 in 2006 (Nõges et al. [68] Figure 6), which is a bit higher than the our estimation for September 2018 (1780 mg C m −2 d −1 , Figure 7a).
There have been large discrepancies in the PP results obtained in different studies using different methods. The long-term study ) by Nõges et al. [68] estimated the average PP from June to October to be 880 mg C m −2 d −1 , based on 14 C methods, while another study [19] estimated the PP from oxygen measurements to be 1632 mg C m −2 d −1 (June to October 2011-2012). Our results would be placed in the middle of those two past studies, with the PP estimated to be 1188 mg C m −2 d −1 from June to October 2018. Our calculations are based on the mean of the entire lake, while the two past studies made their estimations based on measurements from a single point.
PP studies in lakes using remote sensing data have not been very abundant so far. The Great Lakes Primary Productivity Model (inputs are Chl a, K d,490 , q PAR , and also the photosynthesis-irradiance index) to study PP in the Great Lakes (Superior, Huron, and Michigan) using satellite data has been applied in multiple studies [33,74,75]. These already showed promising results of spatiotemporal PP studies and stated the advantages of PP modeling with satellite data. A similar model was applied to remote sensing data by Bergamino et al. [76] to study Lake Tanganyika and Deng et al. [23] to study Lake Taihu. Lake Geneva has been also investigated with remote sensing data using the same model as in the current study [42].
Fahnenstiel et al. [33] estimated the PP year in the Great Lakes to be 5300-8100 Gg C y −1 (Table 2) [75]. Although the sizes of lakes Superior, Huron, or Michigan are not comparable to the smaller lakes targeted in the current study, the PP year,aver showed, for example, that Lake Burtnieks can have more than threefold higher production per m 2 ( Table 2). The lake size plays an important role in the total lake productivity. It is clear that lakes with otherwise similar areal productivity (Lubans, Võrtsjärv, and Tanganyika) have very different PP year , ranging from 18 to 7651 Gg C y −1 ( Table 2). Table 2. The comparison of different lakes by their surface area (S, km 2 ), yearly total primary production (PP year , Gg C y −1 ), and annual average daily areal productivity (PP year,aver , mg C m −2 d −1 ). The oligo-mesotrophic Lake Geneva shows a higher PP year,aver than eutrophic lakes Võrtsjärv and Lubans, but two factors should be kept in mind. Firstly, the PP year,aver values are derived from the PP(z) values over the depth of the euphotic zone, where photosynthesis can occur. In Lake Geneva, the euphotic mixed layer zone can reach up to 30 m (according to Luhtala et al.'s [77] algorithm from the maximum measured SD [78]), while in shallow lakes like Lubans and Võrtsjärv, the euphotic zone cannot exceed the maximum depths of 3.5 m and 6 m, accordingly. Secondly, the lakes are exposed to different climate conditions: Lake Geneva does not have any ice cover. This allows higher production during the entire year, leading to a higher PP year,aver . This might also be one of the reasons, besides the high amount of available nutrients, why Lake Taihu is seemingly the most productive lake, with the highest mean PP year,aver (Table 2). Actually, in the highest productivity month, August, the PP day,aver was 1675.5 ± 710.4 mg C m −2 d −1 [23] in Lake Taihu, whereas in Lake Burtnieks the PP day,aver was 3137.6 ± 1105 mg C m −2 d −1 in August 2018 (Figure 7a).

Lake
Although initially it might seem from the Table 2 that smaller lakes would not contribute significantly to the carbon cycle, lakes with surface areas of 1-1000 km 2 cover about 36% (~1.9 million km 2 ) of the total area of all lakes in the world [5]. Fee et al. [79] stated that the areal production is generally higher for medium-sized lakes, which have a combined offset of temperature and turbulence supporting increased phytoplankton growth, but larger lakes have higher yearly production due to a prolonged growing season. And yet, still a relatively small portion of smaller lakes are monitored regularly or not monitored at all, hence a large portion of the global freshwater productivity is not considered.
The method used in the current work could be easily applicable to larger regions to cover more lakes with different optical properties. The medium-resolution OLCI data where one pixel covers approximately 0.09 km 2 is only suitable for monitoring lakes with a diameter of more than 1 km, leaving the most abundant lake group out of its capability. In addition, although even OLCI data show some spatial patterns, data from high-resolution Sentinel-2 (pixel size from 10 × 10 m to 60 × 60 m, depending on the spectral band) might be even more useful for PP studies of smaller lakes and their spatial variability. This is an especially strong advantage in vegetation-rich lakes like Lubans, allowing the distinguishing of patches of vegetation with higher accuracy, whilst Sentinel-3 pixels are affected by vegetation and thus invalid. Also, Sentinel-2 is proven to have higher quality flagging than Sentinel-3 OLCI [37,72], hence the PP results can be improved. In addition, the PP estimation scheme developed in the current work could be applied to even larger regions and on a longer time scale for interannual comparisons. Even smaller lakes could be included with higher resolution remote sensing data, but in this case, the PP results of the shoreline pixels must be handled with some caution due to the possible influence of the bottom and the adjacency effect; therefore, some additional masks or methods might need to be applied.

Conclusions
The aim of this study was to estimate spatial and temporal PP variability in four lakes by using a PP model based on Sentinel-3 OLCI data. From our studied lakes, only one lake had previous records of PP (Lake Võrtsjärv). Soomets et al. [42] demonstrated that the model gives good results also in deep lakes with a mixed euphotic zone. Consequently, we can also apply the model to the mesotrophic and relatively deep Lake Razna similarly to three other shallow and eutrophic lakes, despite the model being designed for lakes with a well-mixed water column where K d,PAR and Chl a do not change (or only change slightly) with depth [52]. The annual average daily areal productivity, PP year,aver , in 2018 was 333 mg C m −2 d −1 in Lake Razna, 610.2 mg C m −2 d −1 in Lake Lubans, 622.4 mg C m −2 d −1 in Lake Võrtsjärv, and 887.3 mg C m −2 d −1 in Lake Burtnieks. This shows that the most productive lake is the smallest of the studied lakes. The yearly total production of 2018 for the entire lake was 7 Gg C y −1 (Lake Razna), 13 Gg C y −1 (Lake Burtnieks), 18 Gg C y −1 (Lake Lubans), and 61 Gg C y −1 (Lake Võrtsjärv). Even if lake size plays a significant role in the total PP of the lake, the abundance of small and medium-sized lakes would sum up to significant carbon fixation. Hopefully, the free access to Sentinel-2 high-resolution remote sensing data and improved methods will bring the freshwater PP more into the focus of the carbon cycle studies.