Mapping Potential Timing of Ice Algal Blooms From Satellite

As Arctic sea ice and its overlying snow cover thin, more light penetrates into the ice and upper ocean, shifting the phenology of algal growth within the bottom of sea ice, with cascading impacts on higher trophic levels of the Arctic marine ecosystem. While field data or autonomous observatories provide direct measurements of the coupled sea ice‐algal system, they are limited in space and time. Satellite observations of key sea ice variables that control the amount of light penetrating through sea ice offer the possibility to map the under‐ice light field across the entire Arctic basin. This study provides the first satellite‐based estimates of potential sea ice‐associated algal bloom onset dates since the launch of CryoSat‐2 and explores how a changing snowpack may have shifted bloom onset timings over the last four decades.


Introduction
The Arctic is warming nearly four times faster than the rest of the planet (Rantanen et al., 2022), leading to dramatic changes across the land, ocean and atmosphere.One of the most striking changes has been the transformation from an ocean once covered year-round by perennial sea ice to one now dominated by seasonal ice (for e.g., Stroeve & Notz, 2018).As a result, many properties of the ice cover are changing, including its thickness (Lindsay & Schweiger, 2015;Sumata et al., 2023) and the amount of snow that accumulates (Stroeve et al., 2020;Webster et al., 2018).These changes have important implications for the entire Arctic marine ecosystem since its functioning is largely driven by the seasonality of the ice cover and its physical properties.Light is one of the critical drivers of marine primary production, acting as a trigger for sea ice algae and phytoplankton blooms (Ardyna et al., 2020;Leu et al., 2015).Changes in the seasonality of the under-ice light regime also influence the vertical migration of zooplankton (Brierley, 2014;Flores et al., 2023), an active biological carbon pump that accounts for 85-132% of sinking organic carbon (Darnis & Fortier, 2012).While in situ studies have improved understanding of the partitioning of incoming light and its relationship to ice and snow properties (Mundy et al., 2007;Veyssiere et al., 2022), quantifying how the bottom-ice light regime is changing on a pan-Arctic scale over long time-periods requires the use of satellite observations or output from climate models.Stroeve et al. (2021) demonstrated an approach to take existing satellite products to map under-ice Photosynthetically Active Radiation (PAR) from autumn to spring.However, a key limitation was the monthly time-scale, limited by the lack of a daily sea ice thickness (SIT) record.The study was also limited to the cold season due to challenges in processing radar altimetry data during active melt.To handle these constraints, Lim et al. (2022) used weekly ice age data as a proxy for thickness to map the algae habitat period.
Recently Gregory et al. (2021) developed a method to combine radar altimeter data from CryoSat-2 (CS2) and Sentinel-3 (S3), improving the overall temporal resolution with which a daily thickness product can be generated.The same method can be applied to an individual satellite with the caveat that it relies more on temporal interpolation than filling gaps with actual data.Additionally, advances in processing radar freeboards in summer (Landy et al., 2022), and availability of ICESat-2 (IS2) data allow for extension into the melt season.With these new data sets, it is now possible to estimate daily under-ice PAR, and with that, the timing for when algal blooms within and beneath sea ice may occur.
This study expands on our previous work to provide estimates for the potential timing of the annual ice algal bloom onset over the entire Arctic Ocean from 2011 to 2022.We say potential since the study is focused solely on light availability at the bottom layer of the ice where most sympagic algae are concentrated.Though light is not the only requirement, Pinkerton and Hayward (2021) found light transmission explained 69% of the variability in ice algal production in the Southern Ocean.Further, given the importance of snow in limiting light penetration (Perovich, 2007), we additionally assess how the thinning snowpack may have shifted ice algal bloom onset dates since the early 1980s.

Satellite Data Products
We rely on several well-known data sets for the sea ice environment, including passive microwave-derived sea ice concentrations (SICs) provided in near-real-time by the National Snow and Ice Data Center (NSIDC) (Fetterer et al., 2017), and surface albedo, skin temperature and incoming solar irradiance from the Advanced Very High Resolution Radiometer (AVHRR) Polar Pathfinder Extended (APP-X) product (Key et al., 2019).While satellitederived snow depths are available (see Zhou et al., 2021 for a review), we instead rely on modeled estimates from SnowModel-LG (Liston et al., 2020).SnowModel-LG incorporates satellite-derived SIC and ice motion vectors and is forced using ERA5 (Hersbach et al., 2020) or MERRA2 (Gelaro et al., 2017) atmospheric reanalysis.Here we show results using MERRA-2.All data are available daily and at 25-km spatial resolution on the equal-area scalable Earth 2.0 (EASE2) grid (Brodzik et al., 2012).
The next product combines CS2 with Sentinel-3A (launched 16 February 2016) and Sentinel-3B (launched 25 April 2018) (hereafter CS2S3) and processed using both CPOM and LARM waveform retrackers for 2018-2019 and 2019-2020.We are limited to these two years because the L1A waveforms were not available at the time of this study.Our final SIT data set relies on IS2 freeboards from ATL10.
Following Nab et al. (2023), all along-track radar and laser freeboards are first binned onto the 25-km EASE2 grid, with only grid cells containing SIC ≥ 15% selected.Afterward, daily-resolution pan-Arctic radar and laser freeboards are obtained by using a 9-day moving window to fill in data gaps and then converted to ice and snow freeboards, respectively.Finally, SIT is computed assuming hydrostatic equilibrium.More details are provided in the Supporting Information S1.
We considered extending the CS2 record through summer following Landy et al. (2022), but the data were too noisy to produce reliable daily freeboards via our optimal interpolation approach.Instead, we rely on IS2 in summer for 2019 to 2022.We note there may be significant uncertainties in IS2 freeboards during melt, introduced, for instance, by challenges separating melt ponds and leads in the photon classification algorithm, and from variable wind-wave roughening of ponds (Tilling et al., 2018).While absolute magnitudes remain uncertain, our focus is on the relative inter-annual differences leading to variations in bloom timing.
An intercomparison of the daily pan-Arctic SIT for selected winters is provided in Figure S1 in Supporting Information S1.During 2018/2019, CS2 LARM more closely matches IS2 while CPOM SITs are on average 25 cm thicker, with or without the addition of S3.CPOM SIT is also thicker than LARM in 2019/2020 and 2020/2021 yet agrees better with IS2 during freeze-up; later in winter LARM more closely matches IS2.Given that LARM and IS2 agree better during the time when light becomes available, the results below rely on the CS2 LARM data.

Ice Thickness and Snow Depth Distributions
Satellite-derived SIT and modeled snow depth provide mean values at coarse spatial resolution, and thus we must consider sub grid-scale thickness and snow depth distributions for more accurate bottom-ice light estimates.This is particularly important in the marginal ice zone (MIZ), which will have sub grid-scale regions with thin snow and ice and thus considerably larger light penetration.
We keep the ice distribution the same as in Stroeve et al. (2021), but implement an updated representation for snow based on analysis of snow depths from Russian drifting stations and various field campaigns (Mallett et al., 2022).Briefly, a skew normal distribution is used to convert from a mean value into a snow depth distribution using seven discrete bins (Figure S2 in Supporting Information S1).Since a skew normal distribution results in a small fraction of negative snow depths (0.1\% for snow depths of 0-50 cm); negative values are set to 0 cm.This distribution results in overall less light transmitted for deeper snowpacks than the distribution used previously (Figure S3 in Supporting Information S1).

Under-Ice Light Transmission
In modeling light attenuation, several approaches are possible, from exponential decay using empirically derived extinction coefficients for snow and ice (for e.g., Castellani et al., 2017;Lebrun, 2019), to more complex radiative transfer models that rely on inherent optical properties that are directly related to physical state variables (Briegleb & Light, 2007;Stamnes et al., 2018).In the exponential decay approach, the downward solar irradiance is reduced by the surface albedo, and further reduced by an empirically derived extinction coefficient that scales the magnitude of the log-linear decrease of light per unit of thickness.Implementation is relatively straightforward, yet it is not possible to treat multiple scattering or gradients in the inherent optical properties of the snow/ice system.To address the vertical gradient in scattering, exponential approaches can consider including layers such as the surface scattering layer (SSL) that generally covers the top 1-10 cm of the sea ice and results in a large amount of light extinction (Light et al., 2008).
Here we use exponential decay but modify the assumed SSL thickness used previously (10 cm for ice, 3 cm for snow).While there are few direct observations on SSL formation and its seasonal evolution for thin ice and snow, we balance the attenuation in the SSL in proportion to its thickness.For snow depth ≤3 cm, we consider no SSL but use a snow extinction coefficient K s = 40 m 1 .Since observations indicate that the SSL in sea ice needs a drained layer below (Light et al., 2015), we have further introduced a minimum SIT of 50 cm for which an SSL is assumed to exist, and a corresponding linear transition of the SSL thickness from 0 to 10 cm up to a SIT of 80 cm (see also Castellani et al., 2022).The introduction of this progressively increasing SSL thickness increases the amount of PAR transmitted through snow-free thin ice, especially in the MIZ where the ice is often thinner than 80 cm.We also update the extinction coefficients based on Lebrun et al. (2023) and split the snow into two categories depending on if the snow is wet or dry (Table S1 in Supporting Information S1).
While we investigated using the CESM3 1-D Delta Eddington (DE) approximation (see: https://github.com/andypbarrett/seaice_radiative_transfer), the difference between approaches is not large for SIT > 1 m (Figure S4 in Supporting Information S1) and given the order of magnitude time difference in the computation (0.15 ms per loop for exponential decay vs. 1.54 ms for DE) we opted for the exponential decay approach to scale up to the pan-Arctic.In the MIZ, we may be too conservative in our light transmission, because except for albedo reduction, it does not account for the loss of energy outside the asymptotic regime.Another compelling factor bolstering our decision is that the literature contains a richer body of research on extinction coefficients for different sea ice types compared to data pertaining to their physical or inherent optical properties.
Transmitted light (F t , W m 2 ) is converted to PAR (Q PAR , μmol photons m 2 s 1 ), where Q PAR = 0.79 × 4.44F t ) (Stroeve et al., 2021), and PAR transmission is only computed to the bottom of the sea ice; PAR through open water is not included.To scale the satellite albedo which contains open water and sea ice, we use a modified albedo that is weighted by the SIC and using an open water albedo value of 0.07.To determine the date of algal bloom onset, defined as the compensation light intensity that permits photosynthesis to start, a PAR threshold above 1.78 μmol photons m 2 s 1 for four consecutive days is used.This threshold value is taken from Castellani et al. (2017) and is a mean of values identified from 4 different latitudinal bands.Requiring 4 consecutive days accounts for the characteristic lag in algal growth following illumination prior to the exponential phase and aligns with the 3-4 days lag duration documented previously (Forgereau et al., 2023).

Results
Figure 1 shows ice algal bloom onset dates from 2011 through 2022 for: (a) CS2 LARM (2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2021)(2022), (b) blended CS2S3 LARM (2019-2020) and (c) extended with IS2 (2019-2022).Despite CS2 and CS2S3 thickness data being limited through the end of April, there appears to be enough light to initiate a bloom throughout most of the Arctic before May, with only a few locations near the pole with insufficient PAR, the size of which varies from year-to-year (e.g., in 2013 nearly the entire Arctic Ocean has enough light for a bloom to occur before May, whereas in 2017, areas near the pole did not have enough light).
Not surprisingly the more southerly regions experience earlier blooms, with regions outside of the central Arctic Ocean having enough light to initiate ice algal growth already in February to mid-March.Using NSIDC regional masks (see: https://nsidc.org/data/nsidc-0780/versions/1),bloom onset in Baffin Bay occurs on average March 3-5 (Table S2 in Supporting Information S1).Regions within the Arctic Ocean have average bloom onset dates that range between March 11 (Chukchi Sea) and March 27 (Barents Sea); in the central Arctic bloom begins April 3.The addition of IS2 fills in regions that had not yet detected enough light to initiate a bloom, leading to later regionally-averaged bloom onset dates that can be more than a week later in an individual year/region.This stresses the importance of having thickness data that extend into the melt season.
Interannual variability in bloom onset is largely determined by snow depth (Figure S5 in Supporting Information S1).In April 2017, deep snow stretched across the pole and limited light penetration, resulting in large regions with insufficient light by the end of April.Similarly, in 2019, thick snow north of Greenland limited light penetration through April yet extending the time-series with ICESat-2 allowed us to capture a potential bloom onset in May.Incorporation of IS2 also allowed for bloom onset to be detected within the northern areas of the Canadian Arctic Archipelago (CAA).
While the altimeter SIT data record is not long enough to capture trends in bloom onset, given that snow is the dominant factor limiting light penetration (Perovich, 2007), we are interested in how snowpack changes may be driving bottom-ice PAR changes, and thus the sympagic bloom onset.Snow depth has been declining throughout much of the Arctic Ocean since the early 1980s (Figure 2, top panels), with the largest negative trends in the MIZ, increasing in magnitude from January to June.Negative trends are statistically significant in the western Beaufort, East Siberian, Kara and Barents seas.Small positive trends north of Greenland and stretching toward the pole are found but are largely not statistically significant at the 95% confidence interval.Similarly, negative trends in the East Greenland Sea, the Bering, southern Chukchi and eastern Beaufort seas are not statistically significant.On the other hand, statistically significant trends toward increased surface albedo dominate the central Arctic Ocean in March and April, before turning negative in May and June (Figure 2, middle panel).While these albedo trends are small, they would reduce light transmission at the time of year critical for bloom onset.Note, trends in incoming solar irradiation are mostly not statistically significant except the positive trends in March in the central Arctic Ocean and the negative trends in April in the southern portion of the Barents Sea (Figure 2, bottom panel).
To examine how longer term changes in the sea ice environment may influence bloom timing, we re-ran the under-ice PAR calculations starting in 1982 using SnowModel-LG snow depths together with SIT, albedo and incoming solar radiation from APP-X.Wang et al. (2016) previously found APP-X SIT agreed within 48 and 20 cm to those from ICESat and CryoSat-2, respectively.APPX SIT trends can be found in Figure S6 in Supporting Information S1.Two cases are run, one with the Mallett et al. (2022) snow depth distribution and one using the mean snow depth per grid cell.
Over the full time-series, trends toward earlier bloom onset are found throughout large parts of the Arctic Ocean (Figures 3a and 3b).However, in the MIZ where statistically significant reductions in snow depth are found, applying the snow depth distribution does not result in significant changes in the bloom onset (Figure 3a).The reasoning is that in regions with already thinner snowpacks, applying a snow depth distribution allows for "windows" of light penetration already in the 1980s.Instead, in the central Arctic Ocean where deeper snowpacks generally occur, the combination of a thinning snowpack with the snow depth distribution function significantly shifts the bloom onset earlier in the year.Thus, the shape of the snow depth distribution function applied is nontrivial.
To quantify the relative contributions of snow depth, albedo and downwelling irradiance, we compute regionally averaged under-ice PAR using 1982-2010 averages for two of the variables while keeping one varying; SIT and SIC vary as before (Figures 3c-3f).Assuming a mean snow depth per grid cell, we find that snow depth changes play the dominant role in bloom onset trends in the central Arctic Ocean, the Siberian, Laptev and Barents seas.Albedo and downwelling solar irradiance play more important roles outside of the central Arctic and have similar magnitudes within the Arctic Ocean.The notable exception is in the central Arctic Ocean and the Beaufort Sea where downwelling solar irradiance changes largely contribute to the overall changes in bloom timing.

Discussion
To assess how well our approach models PAR is challenging given the lack of in situ data to compare against, especially prior to the melt season.Nevertheless, we evaluate our method against under-ice radiometer data collected during two campaigns to the Arctic Ocean in May-June 2015 and June-July 2017 (Castellani et al., 2020).Using local snow depth, SIT, albedo and downwelling irradiance, we compare computed and measured transmittance (Figure S7 in Supporting Information S1).Results vary by campaign, with a general underestimation in 2017, and more balanced results in 2015.As the 2017 campaign occurred during summer, this underestimation may reflect the lack of implementation of melt ponds in our approach.We additionally evaluated our method against a limited number of radiometer observations mounted to a L-ARM system at 3 ice stations in summer (between 22/06/2020 and 29/07/2020) over first-and second-year ice.Under-ice PAR computed using K i = 1.0 and 1.1 are shown Table S3 in Supporting Information S1, revealing better agreement for a slightly higher ice extinction coefficient.
Studies have suggested the PAR value needed for ice algal blooms to occur can range significantly, with bottomice algal blooms found under comparatively lower (e.g., PAR < 0.17 μmol photons m 2 s 1 , Hancke et al., 2018) or higher (e.g., PAR >20 μmol photons m 2 s 1 , Campbell et al., 2022) light levels depending on algal acclimate state.A lower light threshold may be most realistic given multiple scattering within the sea ice that can increase cellular absorption of light (Ehn & Mundy, 2013) as well as the shade adaptations of ice algae (Johnsen & Hegseth, 1991).While it remains unclear what PAR threshold most broadly represents that needed to initiate ice algal photosynthesis on a pan-Arctic scale, our results suggest a PAR threshold of 1.78 μmol photons m 2 s 1 provides enough light for a bloom to develop already in March across large parts of the Arctic Ocean.Although observations remain sparse, this agrees with previous studies (Campbell et al., 2017;Leu et al., 2015).While we generally cannot resolve bottom-ice PAR in the CAA due to the coarse spatial resolution of the ice motion vectors, in 2021 and 2022 a change in data processing allowed for Lagrangian tracking of ice parcels, and hence SnowModel-LG snow depths to be retrieved.With the addition of IS2 data, this allows for detection of ice algal bloom timings within the northern parts of the CAA in May.In more southerly areas, such as Dease Strait, biological buoys in 2017 and 2019 showed ice algal bloom development on 16 March and 12 March, respectively (Yendamuri pers.comm.), in agreement with our timings.
Bottom-ice PAR however is not the only factor influencing bloom timing.The composition of the initial algal population affects the early stages of growth, as the diversity of ice algae results in a range of light adaptations and growth characteristics (Rozanska et al., 2009).Although the seeding mechanisms of ice algae blooms remain unclear, it has been hypothesized that the algal community originates from the uptake of organisms in seawater during ice growth (Gradinger & Ikävalko, 1998) and that selection processes during winter could play a role in shaping the initial community assemblages of the spring bloom (van Leeuwe et al., 2018).Winter observations suggest ice formation timing impacts algal abundance in new ice, with later freeze-up leading to ice growth from seawater of low phytoplankton concentration, resulting in less accumulation of ice algae (Niemi et al., 2011).As a result, trends toward later freeze-up (Stroeve & Notz, 2018) could cause a decrease in overwintering biomass and early spring algal populations, potentially delaying the onset of the spring bloom.
Incorporation of ICESat-2 during the cold season could improve temporal resolution and reduce the reliance on the optimal interpolation scheme to fill in missing data.However, challenges remain in blending radar and laser freeboards, leading to biases between laser-and radar-derived SIT fields.On the other hand, during the time of year when an ice algal bloom occurs, snow still largely covers the sea ice.Previous sensitivity analysis to ±20% SIT uncertainty gave the lowest PAR errors compared to other variables or parameterizations (Stroeve et al., 2021).Figure S8 in Supporting Information S1 confirms little sensitivity on the bloom onset to our choice of thickness product.Further, using only IS2 data appears to have minimal influence on the overall bloom onset dates, though in the regions where there are differences, the bloom onset generally occurs earlier (Figure S9 in Supporting Information S1).
Instead, more attention is needed on how best to represent snow on sea ice, including its distribution over large grid cells.While SnowModel-LG snow depths compared well with independent observations from buoys, magnaprobes and Operation IceBridge (Stroeve et al., 2020), more dedicated field and airborne campaigns are needed to scale up point measurements to a larger satellite pixel and its dependence on different ice types.The Mallett et al. (2022) distribution function was primarily derived for multiyear ice and its applicability to first-year ice is unknown.Its representativeness for a full 25-km grid cell is also unclear as it was limited by the length of the snow depth transects which were on the order of a kilometer long.Overall, the Mallett et al. (2022) snow depth distribution leads to slightly later bloom onset than the homogeneous distribution used previously (Figure S10 in Supporting Information S1).
We concentrated our efforts on examining large spatial scales and advanced current knowledge using two key considerations: (a) addressing the propagation of light in snow using probability density functions of snow thickness.This approach is crucial for obtaining realistic estimates of under-ice PAR because of the non-linear reduction of light with snow depth.In particular, before the melt season has begun, the density, grain size, depth and thermal history of the snow cover play an outsized role in the amount of light reaching the underside of the ice, and the impact of ice thickness and optical properties remains of second order; (b) leveraging novel breakthroughs in satellite remote sensing observations to infer ice thickness during the crucial period when snow melts and becomes less opaque, leading to a relatively larger impact of ice on under-ice PAR availability.This unprecedented use of satellites benefits our understanding of the changing Arctic Ocean.We stress that joint probability density functions for snow and ice thicknesses evolving over the melt season as well as validation data under large expanses of sea ice would help improve modeling of light transmission.Finally, we note that the APP-X albedo can be biased high for snow-covered sea ice (e.g., values exceeding 0.90), resulting in bottom-ice PAR values likely conservative at times.To reduce these uncertainties, future efforts should incorporate PAR albedo products (Laliberté et al., 2022).
Our work could be further extended to phytoplankton blooms, as ice algal sloughing from the sea ice into the upper ocean can further seed pelagic production (Selz et al., 2018), and the timing of snow melt and pond formation is crucial to assess timings for sympagic aggregate formation (Fernández-Méndez et al., 2014) and phytoplankton bloom evolution (Arrigo et al., 2012;Ardyna et al., 2020).Mapping of melt pond fractions and depth remains problematic from satellite, yet efforts to apply machine learning to melt pond fraction mapping (Lee et al., 2020) from optical data as well as depth from IS2 (Buckley et al., 2023) could pave a path forward to include melt ponds in under-ice PAR calculations.However, while optical imagery covers the entire Arctic daily, Geophysical Research Letters 10.1029/2023GL106486 IS2 requires a full month to cover the entire Arctic, limiting its usefulness for mapping daily pan-Arctic pond depths.
The most prominent challenge remains the heterogeneity of sea ice at smaller spatial scales.The satellite data we used in this study provide large-scale information on key variables needed to model light transmission, yet cannot resolve other important considerations in light attenuation, such as the presence of absorbing particulates like sediment and biological material entrained in the ice (including the presence of algae) or brine/gas inclusions, such that approaches which consider small-scale heterogeneity should be emphasized in future research.

Summary
The Arctic Ocean is experiencing an unprecedented transformation, with implications on the entire marine ecosystem.While much remains unknown on how biodiversity and ecosystem functioning are changing, today's sea ice habitat is characterized by a thinning snow cover, leading to greater light availability for sympagic microorganisms.Satellite data on key physical sea ice properties, together with radiative transfer modeling provides an approach to evaluate how the phenology of bottom-ice light availability is changing, and as a result, ice-algal blooms, a critical food supply for grazing organisms (Kohlbach et al., 2016).Light has other implications as well, it affects the food quality of ice algae, such that under more light, sea ice diatoms produce less lipids (Leu et al., 2016).Furthermore, a thinning snow and sea ice cover also have important implications for timing of ice algal bloom termination (Campbell et al., 2015) that require further assessment.Finally, the true development of a bloom implies sustained growth over a period of time that requires not only light, but other environmental conditions to be met.This study provides a method by which daily under-ice PAR could be combined with other information to assess changing primary productivity in the Arctic as the sea ice cover and its overlying snow cover continue to thin and decline.

Figure 1 .
Figure 1.Yearly potential ice algal bloom onset dates from 2011 to 2022.Blended CryoSat-2 and Sentinel-3 data are for 2019 and 2020 only, extensions with ICESat-2 are outlined in red.Areas in white near the pole have insufficient light to initiate a bloom.In 2021 and 2022, the APP-X data set has a circle of missing data around the pole.

Figure 2 .
Figure 2. Monthly averaged snow depth (top), surface albedo (middle) and Incoming Solar Irradiance (bottom) trends from January through June between 1982 and 2018.Hatches denote trends statistically significant at the 95% confidence.

Figure 3 .
Figure 3. Trends in bloom onset-date dates from 1982 to 2018 using (a) the Mallett et al. (2022) snow depth distribution and (b) mean snow depth; (c)-(f) contribution to bloom onset trends due to varying snow depth, albedo, incoming solar irradiance and all varying respectively.Results for (c)-(f) are shown as regional averages for each National Snow and Ice Data Center region.