A Multi-Sensor and Modeling Approach for Mapping Light Under Sea Ice During the Ice-Growth Season

Arctic sea ice is shifting from a year-round to a seasonal sea ice cover. This substantial transformation, via a reduction in Arctic sea ice extent and a thinning of its thickness, influences the amount of light entering the upper ocean. This in turn impacts under-ice algal growth and associated ecosystem dynamics. Field campaigns have provided valuable insights as to how snow and ice properties impact light penetration at fixed locations in the Arctic, but to understand the spatial variability in the under-ice light field there is a need to scale up to the pan-Arctic level. Combining information from satellites with state-of-the-art parameterizations is one means to achieve this. This study combines satellite and modeled data products to map under-ice light on a monthly time-scale from 2011 through 2018. Key limitations pertain to the availability of satellite-derived sea ice thickness, which for radar altimetry, is only available during the sea ice growth season. We clearly show that year-to-year variability in snow depth, along with the fraction of thin ice, plays a key role in how much light enters the Arctic Ocean. This is particularly significant in April, which in some regions, coincides with the beginning of the under-ice algal bloom, whereas we find that ice thickness is the main driver of under-ice light availability at the end of the melt season in October. The extension to the melt season due to a warmer Arctic means that snow accumulation has reduced, which is leading to positive trends in light transmission through snow. This, combined with a thinner ice cover, should lead to increased under-ice PAR also in the summer months.


INTRODUCTION
The Arctic is undergoing a period of profound transformation in response to anthropogenic warming, with the loss of the sea ice cover one of its starkest changes. Since the late 1970s, the summer ice cover has shrunk in area by 40%, whereas changes in winter have been much smaller, on the order of 10%. This loss in sea ice area has been accompanied by loss of the thick multiyear sea ice, which today makes up only 30% of the Arctic Ocean sea ice compared to 70% 40 years ago (Maslanik et al., 2011;Stroeve and Notz, 2018). These changes reflect thinning of the ice cover, which has thinned over the central Arctic basin by 65% since 1975 (Lindsay and Schweiger, 2015). While the 2012 sea ice extent minimum was striking, variability is increasing, in relation to a thinning ice cover (Goosse et al., 2009) and recent departures from average conditions during the transition seasons have become even more anomalous than in summer: Arctic sea ice extent (defined as the area with at least 15% sea ice concentration, SIC) in May and November 2016 fell nearly 4 standard deviations below the 1981-2010 long-term average (Stroeve and Notz, 2018). This represents a distinct change in the seasonality of the Arctic Ocean as the melt onset is happening earlier and the fall freeze-up later (Stroeve et al., 2014;Stroeve and Notz, 2018;Lebrun, 2019).
Shortwave radiation, and its visible fraction in particular, provides an essential control on the Arctic surface energy budget (Maykut, 1986;Perovich et al., 2007) and on microbial ecosystems (Wassmann and Reigstad, 2011). In addition, light transmitted under sea ice warms the upper ocean and in turn drives basal sea ice melt (e.g., Maykut and McPhee, 1995;Vivier et al., 2016).
In regards to ecosystems, light is the main energy source for the development of phytoplankton and sea ice algae. Together sea ice algae and phytoplankton form the base of the Arctic marine food web, sustaining directly sea ice associated macrofaunal and pelagic zooplankton (Kohlbach et al., 2016(Kohlbach et al., , 2017. The sea ice scape greatly affects the amount of light reaching the upper ocean. Sea ice and especially its snow cover are excellent reflectors, whereas melt ponds and open water within the pack effectively act as windows, increasing light supply to the surface ocean (Frey et al., 2011;Assmy et al., 2017). In turn, sea ice drastically attenuates light reaching the ocean surface, which would reduce the depth over which waters are biologically productive (Sverdrup, 1953;Horvat et al., 2017). On the other hand, the ice underside provides a highly variable and heterogeneous habitat for different ice-associated macrofaunal as well as certain zooplankton whose vertical migration is often triggered by food availability and periodic changes in light availability (e.g., Berge et al., 2014).
As the winter ends and the sun returns above the horizon again, light is the major factor controlling algal growth onset (Castellani et al., 2017) and production (Horner and Schrader, 1982). Ice algae are able to take advantage of low under-ice photosynthetic active radiation (PAR) levels during this time of year (Mock and Gradinger, 1999) at sites with a relatively low snow depth cover. Snow distribution in spring (March-April-May) is thus the major physical driver of ice algae phenology, and changes in sea ice and snow scape are expected to affect both phytoplankton and ice algal activity. The clearest example is provided by earlier ice retreat and delays in freeze-up, which have cascading impacts on light availability. For one, earlier melt onset allows light to enter the Arctic Ocean earlier and closer to the summer solstice than it used to, whereas delays in freeze-up have extended the light season into early winter, leading to visible Arctic planktonic activity earlier in spring and later into fall (Arrigo and van Dijken, 2011;Ardyna et al., 2014). The extended open water season not only increases the coupling between the ocean and the atmosphere, but additionally reduces the amount of time over which snow can accumulate on the sea ice (Webster et al., 2014;Stroeve et al., 2020). Other expected changes in sea ice scape, such as melt pond fraction and depth, currently challenging to detect (Zhang et al., 2018), would also influence marine autotrophs in polar seas in summer. Together these changes have important implications for the in-ice and underice biota, influencing light availability, ocean properties, and the timing of sea ice algae and phytoplankton blooms (Bluhm et al., 2017). These changes in sea ice scape are also accompanied by stratification and circulation changes, which also affect under-ice plankton activity, via the modulation of nutrient fluxes, critical to the total possible photosynthesis over the year (Vancoppenolle et al., 2013;Randelhoff et al., 2020). In other words, changes in the sea ice and snow characteristics alters phenology of primary productivity in the Arctic Ocean.
However, our understanding of how this sea ice changes impact primary productivity is still in its infancy. Our understanding of ecosystem function, sea ice, and upper ocean processes in the Arctic Ocean has been mostly derived from a multiyear ice setting, rather than the thinner first-year ice dominated Arctic of recent years. As a result, our current climate models use formulations of light transmission through sea ice in large part based on parameterizations for multiyear sea ice (e.g., Grenfell and Maykut, 1977;Fichefet and Morales Maqueda, 1997;Perovich et al., 2002;Briegleb and Light, 2007). Recent observations have shown that the transition from a multiyear to first-year ice-dominated Arctic Ocean has increased the amount of light reaching the upper ocean, with a threefold increase in light transmittance (Nicolaus et al., 2012). This was primarily a result from increased melt pond fraction over first-year sea ice compared to multiyear ice. Nicolaus et al. (2012) also found that energy absorption in first-year ice was 50% larger.
Quantifying the availability of light under the ice is key if we are to better understand how primary productivity functions in today's Arctic, and how it may change in the future. Specifically, we need to better quantify the availability of PAR (400-700 nm) under the ice, and how this affects ecosystem function. This requires improved parameterizations of the light climate under an Arctic Ocean dominated by firstyear sea ice with reduced winter snow cover and increased melt pond fraction in summer. The most comprehensive set of measurements from six years of in situ observations from spring to autumn were recently published by Katlein et al. (2019). These data, while still somewhat limited in both space and time, are improving our ability to derive robust relationships between ice conditions (e.g., thickness, melt pond fraction, snow depth) and light transmission.
To scale up to the pan-Arctic level, we must rely on satellite observations and model-produced data sets of the key system controls. Satellites provide several sea ice variables needed to quantify under-ice light, including (a) surface albedo, (b) melt pond fraction, (c) timing of melt onset, (d) SIC, (e) sea ice thickness (SIT), and (f) surface topography. A key limitation has been the snow depth on sea ice, which has generally been poorly mapped with satellites. This is particularly important before the melt season starts because snow depth plays a larger role than ice thickness in limiting light availability under the ice (e.g., Mundy et al., 2005). Nevertheless, several recent advances in modeling snow accumulation over sea ice have been made, using a combination of modeling with atmospheric reanalyses and satellite drift datasets (e.g., Blanchard-Wrigglesworth et al., 2018;Petty et al., 2018;Liston et al., 2020), as well as combining dual radar frequencies from satellite altimeters (e.g., Lawrence et al., 2018).
Another key limitation is the lack of ice thickness information in summer, as well as the temporal and spatial coverage during winter. Prior to the launch of ICESat-2, pan-Arctic SIT has mostly been monitored using satellite radar altimeter missions that span 1993 through present. These have been obtained at various temporal and spatial resolutions, but are entirely limited to the cold (non-melt) season (i.e., October through April). Since 2010, CryoSat-2 has provided monthly pan-Arctic SIT estimates (e.g., Laxon et al., 2013;Tilling et al., 2018). With the launch of ICESat-2 in September 2018, SIT may potentially also be retrieved in summer. Further, the combination of ICESat-2 and CryoSat-2 freeboards may give a direct estimate of snow depth (Kwok and Markus, 2018), but this strongly depends on the location of the dominant scattering horizon of CryoSat-2, which is unlikely to be from the snow/ice interface (e.g., Willatt et al., 2011). Other important sea ice features for quantifying light availability under the ice include (1) melt ponds, (2) ridges, and (3) leads. Each of these have been mapped using satellites with either limited success in the case of melt ponds (Tschudi et al., 2008;Rösel and Kaleschke, 2012;Scharien et al., 2017;Yackel et al., 2018), or not mapped over large spatial and long temporal scales in the case of leads (e.g., Eicken et al., 2006;Zakharova et al., 2015;Willmes and Heinemann, 2016;Lewis and Hutchings, 2019), and while surface roughness has been detected (e.g., Landy et al., 2015;Petty et al., 2016;Nolin and Mar, 2019), ridge height has yet to be retrieved from satellites. Nevertheless, ridged ice can be important biological hotspots (e.g., Fernández-Méndez et al., 2018).
Given that current uncertainties in the light limitation level (Popova et al., 2010), as well as nutrient uncertainties undermine confidence in primary production for the 21st century (Vancoppenolle et al., 2013), progress is needed on all fronts. This paper focuses on our capabilities for mapping light between 400 and 700 nm under ice with currently available satellitederived products and state-of-the art light parameterizations, and suggests ways forward to fill important gaps to improve our understanding on how light transmission is changing in a changing sea ice environment. Given the limitations with today's satellite products on a pan-Arctic scale, we are presently restricted to monthly means over the October to April time-period. Because of these inherent difficulties, we focus here on the CryoSat-2 time period (October to April, 2010 to 2018), and examine how accurately we can estimate under-ice light fields in a relevant manner for physical and ecosystem studies at the pan-Arctic scale. Prior to the melt season, snow depth will play the dominant role in limiting light entering the upper ocean (e.g., Katlein et al., 2019), and thus the focus is largely on how snow variability and change impact under-ice light levels. We further extend our timeperiod to May and July to examine how snow depth trends impact light transmission through snow as the snowpack begins to melt.

Solar Transmission Calculations
To map the under-ice light with the satellite products, we rely on a generalization of the approach proposed by Maykut and Untersteiner (1971) and Grenfell and Maykut (1977), based on a specification of apparent optical properties (AOPs), namely surface albedo (α) and vertical attenuation coefficients (κ), assuming a two-level Beer-Lambert exponential decay in the snow-ice system. This approach considers that the snow-ice system is comprised of two layers, which absorb most of the radiation, the so-called single scattering layer (SSL), and a lower layer, where, to a first approximation, the attenuation of light follows the Beer-Lambert Law (e.g., Maykut and Untersteiner, 1971). While the details of the formation of the SSL are not well known, it is observed to persist widely across melting sea ice (Untersteiner, 1961;Perovich et al., 2002;Light et al., 2015) and is thus assumed to form due to surface melting, with a thickness in the range of 1-10 cm (Light et al., 2008). Alternative approaches have also used three layers, by separating the sea ice column below the SSL into a 'drained layer' (DL), and one that is under the waterline (Light et al., 2015, see also section Appropriateness of Methodology). Below, we model the snow-ice system as consisting of a series of 3 (2) layers, if snow is present (absent). On top of the system, we assume a thin, optically defined near-surface highly scattering layer (SSL) of thickness h o . Below the SSL, the remaining energy is absorbed according to the Beer-Lambert Law in the remainder of the snow (if any) and ice. The absorption of radiation in the SSL is described by the surface transmission parameter (i o ) (Maykut and Untersteiner, 1971;Grenfell and Maykut, 1977;Light et al., 2008).
The original parameterization of Maykut and Untersteiner (1971) considers a single surface type. Here we expand this formulation to incorporate various surface types in a single satellite pixel. To accomplish this, we express the transmitted, broadband irradiance (F t , W m −2 ) available at the ocean-sea ice interface, as a function of the downwelling solar irradiance at the atmosphere-ice interface (F o , W m −2 ), assuming a mix of n sea ice/snow/water types over the satellite pixel as a sum of the area fractions (A j ) of each ice type multiplied by the modelcomputed transmittance (T j ) (ratio of transmitted over incoming irradiance): For open water (e.g., leads), the transmittance is T w = (1-α w ), where the albedo of open water (α w ) is typically around 0.07 (Pegau and Paulson, 2001). For other sea ice types, we use the two-level Beer-Lambert approach. In this context, for snowcovered sea ice, transmittance (T s ) reads: where h s is the snow depth and h i is the ice thickness. Similarly, for bare ice the transmittance (T i ) is: Finally, through melt ponds, we get: where α is the albedo, k is the attenuation coefficient, and i o and h o correspond to the transmission parameter and thickness of the surface scattering layer that forms on top of snow, if present, or ice (if no snow is present). Another difference from the original approach of Maykut and Untersteiner (1971) is that we use coefficients that were optimized based on observation-based transmittance data, over the entire shortwave range (0.2-4.0 µm). More precisely, the coefficients used in equations 1-4 were adjusted (Lebrun, 2019) using the GreenEdge dataset Randelhoff et al., 2019;Massicotte et al., 2019a). The GreenEdge field activities took place in 2015 and 2016 within a seasonal ice environment, and were focused on plankton and marine biogeochemistry in Baffin Bay. Two land-fast sea ice camps were set from April to July 2015/2016 near Qikiqtarjuaq Island on the west coast of Baffin Bay. A cruise onboard the CCGS Amundsen across the pack ice edge occurred in June-July 2016 in Baffin Bay. Two hundred and fifty eight combined observations of fast and pack ice thickness, snow depth, melt pond coverage, and above-and under-ice spectral irradiance, were retained and processed to adjust the parameter values summarized in Table 1. Adjustment was made in order to achieve a reasonable match between observed and retrieved shortwave transmittance (Figure 1), keeping the parameters within the uncertainty range, based on available observational studies (Perovich and Gow, 1996;Light et al., 2008;Järvinen and Leppäranta, 2011). As all Green Edge observations were made in seasonal ice, the optimized parameters have to be seen as representative of a first-year sea ice environment. Overall, the retrieved and transmittance show a reasonable match, but uncertainties remain large, especially at low transmittance. Processing and analysis are detailed in Lebrun (2019).
Note, for snow there are two values used, one for a dry snowpack and one for a wet snowpack (i.e., melt has begun). Note also that the thickness of the near-surface high scattering layer is assumed as 3 cm for snow and 10 cm for ice. However, the thickness of snow and ice can be less than these values. i o is the surface transmission, h o is the thickness of the medium and k is the extinction coefficient. . Observed transmittance is the ratio of irradiance under sea ice to that above sea ice, both from radiometers deployed in situ. Calculated transmittance derives from equations 1-4, using observed values of snow and ice thicknesses, air temperature and melt pond coverage, and with input parameters from Table 1.
Observational evidence for attenuation in very thin snow and ice is sparse. To regularize attenuation in very thin snow and ice we assume it to be equivalent to the attenuation in the fully developed near-surface high scattering layer (h 0 = 0.03) but reduced in proportion to the actual thickness.

Conversion From Solar Energy to Visible Quanta
Photosynthesis estimates rely upon photon counts in the visible wavelength range (Q PAR , µE/m 2 /s, 0.4-0.7 µm) rather than upon solar energy in the solar waveband (F SW , W/m 2 , 0.2-4 µm). Both are tightly (practically linearly) related, hence the conversion from one to another is straightforward and often used in Earth System model applications. The relationship between Q PAR and F SW can be decomposed as follows: (1) the conversion from solar to visible energy, (2) the conversion from energy to quanta. Formally, one gets: where F PAR /F SW and Q PAR /F PAR can practically be considered as constants. Visible energy, visible quanta and shortwave energy fluxes are defined as (Morel and Smith, 1974): In temperate marine environments, these relationships are well constrained. The ratio of PAR to solar energy ranges overF PAR /F SW = 0.45-0.50, according to simulations with a radiative transfer model (Frouin and Pinker, 1995;Frouin and Murakami, 2007). A classical value for the visible quanta-toenergy ratio is Q PAR /F PAR = 4.60 ± 0.03 J/µE (right under sea surface) decreasing to an average 4.15 J/µE within interior waters, based on a compilation of observations in temperate oceans (Morel and Smith, 1974). Therefore, over the ice-free fraction of the pixels, we use: In ice-covered environments, however, differences arise because of the presence of sea ice. Indeed, both ratios depend on the light spectral distribution, which is altered by snow and biogenic particles in sea ice. GreenEdge data offers a means to evaluate these relationships in ice-covered environments (Figure 2). GreenEdge spectral irradiance (W/m 2 /nm) observations feature 19 channels, with 10 nm bandwidth. Raw data were interpolated on 1-nm-wide bands, covering a reasonable part of the full shortwave range (320-800 nm). Under fast and pack ice environments, both relationships (Q PAR /F PAR , F PAR /F SW ) are nearly perfectly linear (R 2 = 0.999, 0.994). The quantum-energetic ratio is slightly smaller than for ice-free waters: Q PAR /F PAR = 4.44 ± 0.005µE/J. This is because the spectrum is shifted toward blue / green wavelengths, where photons are more energetic, hence one needs less photons per unit energy. A much clearer difference is that the fraction of energy in the visible is F PAR /F SW = 0.79 ± 0.003, i.e., much larger than in ice-free waters, because of the strong attenuation power of ice and snow for infrared radiation. Based on these considerations, over the ice-covered fraction of the pixels (regardless of the presence of snow or melt ponds), we use:

Data Products
While several satellite-derived sea ice data sets now exist, many do not cover a long time-period, or cover the pan-Arctic region.
In this paper we attempt to map light under the ice with currently available state-of-the art light parameterizations and pan-Arctic sea ice and snow products. Given the lack of pan-Arctic SIT information in summer, an important caveat is that we remain limited to the cold season and monthly means, while albedo, snow depth and ice concentration are all available at higher temporal resolution: daily to twice-daily time-steps. Further, pan-Arctic ice thickness estimates are currently limited to monthly means from satellites such as CryoSat-2. Finally, while some limited estimates of lead fractions and surface roughness exist, they are not available at the same temporal and spatial resolution needed for this study. Thus, in our snow depth distributions discussed below, we assume snow accumulated over level sea ice, though snow is redistributed and sublimated with winds (see Liston et al., 2020). The key data sets used are described below and will provide information on snow depth, SIC, SIT, surface albedo, incoming solar radiation and whether or not the snowpack is wet or dry.

Albedo and Solar Radiation
Visible satellites have provided observations of visible reflectance for several decades. A key sensor, the Advanced Very High Resolution Radiometer (AVHRR) has flown on several NOAA POES satellites, providing the potential to map surface albedo as far back as October 1978. AVHRR carries one visible channel, one near-infrared channel, one mid-infrared and two thermal infrared channels, from which surface albedo and surface temperature can be determined. The AVHRR Polar Pathfinder Extended (APP-X) data set (Key et al., 2019), produces twice-daily (ascending and descending passes) estimates of broadband shortwave surface albedo for both polar regions. Here we use the data corresponding to local time of 14:00 GMT (ascending pass). While the surface can only be viewed during clear-sky conditions, a cloud masking algorithm and radiative transfer modeling is used to provide an all-sky (e.g., clear or cloudy) albedo at every satellite pixel. Further, this radiative transfer modeling is used to infer the daily averaged incoming solar radiation at the ground. Early evaluation of the data product was performed through comparisons with in situ observations collected during the yearlong Surface Heat Budget of the Arctic Ocean (SHEBA) experiment (Moritz et al., 1993). Results were summarized by Key et al. (1997). In short, biases with incoming solar radiation were on the order of 9.8 Wm −2 , whereas surface albedo comparisons revealed a bias of 0.028 for the all-sky albedo. Clear-sky albedo conditions had smaller mean errors. While these errors were previously assessed, it is likely that the bias in albedo is larger than what was reported previously. Data for both variables are provided daily on a 25-km equal-area (EASE) grid.

Sea Ice Concentration and Timing of Melt Onset
Since the late 1970s, there have been a series of multi-frequency passive microwave sensors launched that today provide more than 40 years of brightness temperatures that can be used to map sea ice given the large dielectric contrast between open water and ice. Starting with the launch of the Nimbus-7 SMMR instrument in October 1978, and followed on with successive DMSP SSMI instruments (since July 1987) there have been daily (twice daily for SMMR) observations of the polar regions at a nominal spatial resolution of 25 × 25 km 2 . While several algorithms exist for retrieving SIC from these brightness temperatures (e.g., see Ivanova et al., 2015), we rely on the NASA Team sea ice algorithm . This algorithm is produced in near-real-time by the National Snow and Ice Data Center (Fetterer et al., 2017). Observational values (+) are split between the two landfast sea ice camps in Qiqitarjuak (Qiq'15 and Qiq'16) and the drift ice cruise onboard the R/V Amundsen in Baffin Bay (Amu'16). Linear fits use coefficients corresponding to observation-based values (linear regression coefficient and mean ratio) and to classically used values, corresponding to open water situation (Morel and Smith, 1974, MS74;and Frouin and Pinker, 1995, FP95).
While the accuracy of the SICs is generally high during the cold season, during summer its precision is downgraded due to liquid melt water on the snow and/or ice surface that causes the NASA Team algorithm to underestimate the true SIC as these areas will be interpreted as open water. As we are currently limited to the cold season because of the lack of SIT information once melt begins, the choice of algorithm is less important. Importantly, the sensitivity of the microwave emissivity to liquid water in the snowpack provides a means to flag if the melt has begun using the approach of Markus et al. (2009). Both the melt onset and SIC data sets are provided on a 25-km Polar Stereographic grid. All data sets are re-gridded to a 25-km equal area (EASE) grid.

Snow Depth
Snow depth has not been routinely observed from satellites. Some studies have relied upon passive microwave brightness temperatures used in the SIC algorithms to detect snow depth over first-year ice regions (e.g., Markus and Cavalieri, 1998;Comiso et al., 2003), yet results yield unrealistic snow depths in regions when snow melt has started and along the marginal ice zone . Other efforts have involved using atmospheric reanalysis data to model snow accumulation over sea ice, using simple (e.g., Kwok and Cunningham, 2008;Petty et al., 2018) to more sophisticated snow models (Liston et al., 2018. These have been implemented in either a Eulerian (Petty et al., 2018) or Langrangian framework (Kwok and Cunningham, 2008;Liston et al., 2020;Stroeve et al., 2020).
In this study we rely on SnowModel-LG  to generate snow-depth. SnowModel-LG is a spatially distributed snow-evolution modeling system that has previously been demonstrated to be capable of simulating high-resolution snow accumulation around ridges and snow dunes (Liston et al., 2018), and has recently been applied to map snow depth and density on a pan-Arctic scale using satellite-derived ice motion vectors in a Lagrangian framework Stroeve et al., 2020). In short, SnowModel-LG consists of 4 sub-models: (1) EnBal (Liston, 1995;Liston et al., 1999) calculates surface energy exchanges and snowmelt; (2) SnowPack-ML (Liston and Hall, 1995;Liston and Mernild, 2012) is a multilayer snowpack model that simulates snow depth, density, grain size, and habit (e.g., wind slab, depth hoar), and temperature evolution of each storm-related layer. Density evolves as a function of air temperature, overburden pressure, time, and the presence of blowing snow (Liston and Hall, 1995;Liston et al., 2007). The latest version adds vapor-flux related metamorphism responsible for growing faceted crystals and depth hoar ; (2) SnowTran-3D (Liston and Sturm, 1998;Liston et al., 2007) simulates formation and distribution of snowdrifts that develop around topographic obstructions to wind. It simulates the horizontal snow-transport flux at each time step, as a function of surface shear stress from the wind, surface shear strength of the snow, and amount of snow available for transport; and (3) SnowDunes (Liston et al., 2018) simulates high-resolution snow surface features such as dunes and sastrugi. Ice motion vectors derived from satellite play a key role in redistributing the snow for the pan-Arctic simulations.
In this study we use the daily snow depth and density product generated from MERRA-2 atmospheric reanalysis (Gelaro et al., 2017) together with ice motion vectors from Version 4 of the Weekly EASE grid sea ice motion vector data set (Tschudi et al., 2019). Comparison against various in situ data sets show SnowModel-LG simulations capture the observed spatial and temporal variability of snow accumulation . The data set spans August 1980 through July 2018 and is offered, as with the other products used, on the 25-km EASE grid.

Sea Ice Thickness
While we have more than 40 years of SIC observations from satellite, we do not have a similarly long-term SIT data set. Several satellite missions have provided insight into SIT using either radar (e.g., ERS1/2, Envisat, CryoSat-2) or laser altimetry (ICESat-1/2), though they are not consistent in time or spatial area covered. In this study, we focus on the 9 years of continuous ice thickness observations from CryoSat-2. As before, several different algorithms exist to retrieve SIT from CryoSat-2 (e.g., Kurtz et al., 2014;Tian-Kunze et al., 2014;Hendricks et al., 2016;Lee et al., 2016;Kurtz and Harbeck, 2017;Ricker et al., 2017;Tilling et al., 2018). Note that CryoSat-2 does not measure SIT directly but measures the ice freeboard, which together with snow depth, snow density and ice density can be converted into SIT. In this study we chose to use a new algorithm developed by Landy et al. (2020), which uses a Lognormal Altimetry Retracker Model (LARM) applied to CryoSat-2 returns, together with snow depth and snow density from SnowModel-LG to estimate SIT. The LARM algorithm is based on simulations of CryoSat-2 waveforms performed with a physical model for the SAR altimeter echo backscattered from sea ice (Landy et al., 2019). The physical echo model accounts for realistic variations in the sea ice surface roughness and radar backscattering properties, at the scale of the CryoSat-2 footprint, which affect the derived sea ice freeboard. Since SnowModel-LG has been run with both MERRA-2 and ERA-5 atmospheric reanalysis, an average of both is taken to represent the snow depth and density in the SIT-retrieval. Data are provided monthly on the 25 km EASE grid.

Use of a Heterogeneous vs. Homogeneous Snow Depth and Ice Thickness Distribution
An important issue to be aware of when using mean snow depth and ice thickness products for light transmission is that at 25km grid scales, observations clearly show that snow depth and SIT do not follow simply mean values per grid cell (i.e., h i and h s ). They instead follow a probability density distribution (PDF) with varying widths of distributions depending on time of year and ice type (e.g., Renner et al., 2013). Since thin ice or thin snow cover play an important role in light transmission we would underestimate light transmission by simply assuming mean values for h s and h i . In terms of SIT, models often parameterize the sub-grid-scale SIT variations by replacing h i by an ice thickness distribution (Thorndike et al., 1975) that has been discretized using fixed PDFs per thickness category (Hibler, 1984). However, several studies have highlighted the fact that ice thickness may have more than one mode, with bimodal distributions found in regions with both first-year and multiyear ice present (e.g., Haas et al., 2010). Based on airborne electromagnetic induction sounding (EM-bird) measurements reported in Haas et al. (2010), Castro-Morales et al. (2014) evaluated the impact of using 15 ice categories on the surface heat budget in the Arctic Ocean. This representation of ice thickness resolves more thin ice than the flat (e.g., same probability) seven thickness categories used in some sea ice models. Because not enough is known as to how ice thickness distributions may be changing as the ice cover has thinned and become more mobile (e.g., Rampal et al., 2009), we decided to model the ice thickness distribution through 15 ice thickness categories following Castro-Morales et al. (2014), such that the ice thickness is distributed between 0 and a thickness of 3h i with a bin width of 3h i /15 = 0.2h i .
As we have for SIT, snow depth is highly heterogeneous at small spatial scales. Abraham et al. (2015) show for example that, for the same mean snow depth, this small-scale heterogeneity tends to increase transmission below sea ice because of the large contribution of thin snow. They tested various distributions (gamma, Rayleigh) and compared them to the assumption of a snow distribution between 0 and twice the mean snow depth, showing little difference. Thus, in this study we also assume a snow distribution similar to ice thickness but instead solve for the integral of the probability density function g(h s ). For example, if we just consider the snow transmittance and assume snow is homogenous with mean thickness, the snow transmittance can be given by: To account for the fact that snow is heterogeneous, characterized by g(h s ) we can write: where A is snow concentration or area (fraction of grid-cell covered by snow). The first moment of g(h s ) is the mean thickness: Using, the first integral above (Eq. 11) one then gets: Calculating transmittance with g(h s ), we get: Many models assume that thin ice has a thinner snow pack, yet we lack the observations to fully support this assumption.
Since we do not fully understand how snow depth distribution depends on ice thickness distribution we instead solve the full snow distribution per ice thickness category. Thus, including the ice thickness per thickness category h i,n , we obtain: We conclude with a discussion on uncertainty of satellite and transmission parameters on the results. Specifically, for each month we additionally estimate locally the error in under-ice light transmittance as a function of uncertainty in the physical sea ice/snow parameters (ice thickness, snow depth, surface albedo, and SIC) as well as the attenuation properties ( Table 2). These uncertainties are considered with respect to errors in the satellite retrievals (snow depth, albedo, ice thickness, and ice concentration) as well as the parameters used in the light calculations and are based on current knowledge. For the albedo, we assume an uncertainty of 0.05, higher than the accuracy reported in the APP-X data set, as a result of problems in handling the anisotropic reflectance of snow, particularly under large solar zenith angles (i.e., February). At the ice edge however, the errors may be larger given the fact that the ice edge is not well-resolved from passive microwave satellite data, and thin ice limitations from CryoSat-2 become important. This is discussed further in section "Appropriateness of Methodology."

Pan-Arctic Under-Ice Light
Monthly mean (February to April and October) summaries of under-ice PAR, transmittance and the various input data from 2011 to 2018 are summarized in Figures 3-9, respectively. All of these results include both the snow and ice distributions as described in Section "Data Products." The impact of using our modeled snow depth and thickness distributions compared to mean values per grid cell are discussed in section "Impact of Homogeneous vs. Heterogeneous Snow Depth and Ice Thickness Distributions." Some key results stem from these maps. First, in April there is a substantial amount of under-ice PAR (>10 µmol m −2 s −1 ) entering the Arctic ocean in the southerly locations (Figure 3).  Grenfell and Maykut, 1977. This is expected as there is more incoming solar radiation in these locations compared to the central Arctic and seasonally more in April than the other months evaluated. There is also less snow in these regions, often less than 8 cm, which combined with incoming solar radiation in excess of 400 Wm −2 , allows for a larger light transmission (Figure 4) and therefore, under-ice PAR. Second, interannual variations in the under-ice PAR are strongly linked to the snow depth variability (Figure 5). This is particularly true from February to April, whereas SIT largely contributes to October under-ice PAR variability (Figure 6). The lowest pan-Arctic snow accumulation in April for example occurred in 2013, resulting in an increase in under-ice PAR over larger parts of the central Arctic Ocean compared to other years (e.g., 4-10 µmol m −2 s −1 compared to on average 2-5 µmol m −2 s −1 ). Conversely, in 2011, snow depths across large parts of the Arctic Ocean exceeded 35 cm, which in turn limited the underice PAR to values less than 10 µmol m −2 s −1 most everywhere in the central Arctic. Spatially, interannual variability of underice PAR is large in the Beaufort Sea, reflecting large variations in snow depth in this region (e.g., compare 2016 to 2011).
While there is less light in February and March, thin snow and ice thickness in regions such as the Bering Sea, Sea of Okhotsk and Baffin Bay already allow for increased light transmittance and under-ice PAR values exceeding 20 µmol m −2 s −1 . Similarly, in some years, under-ice PAR in the Chukchi and Beaufort seas are seen to increase to above 10 µmol m −2 s −1 already in February. Hancke et al. (2018) have recently showed that net growth in ice algae could be initiated at irradiances lower than 0.17 µmol m −2 s −1 . Based on other studies, threshold values for onset of algal growth are in the range of 0.5-5 µmol m −2 s −1 (see Letelier et al., 2004;Tremblay and Gagnon, 2009;Leu et al., 2015). These maps suggest that light levels are higher than the lower threshold value in February for almost all regions of the Arctic Basin. This suggests that algal growth can start even at this early state in the season. If one considers a higher threshold than 0.17 µmol m −2 s −1 (e.g., Hancke et al., 2018), we find that, algal growth would happen between February and March at lower latitudes, then move northward in April, without reaching regions further north than 85 • N. Only in 2012 and 2013 did light levels in the     high Arctic (above 85 • N) become large enough to initiate algal growth in April. How do these results compare to published values? In the months of February to April, transmittance values remain below 0.01 in most of the Arctic, except for the shelf areas, where thinner ice and snow cover allow for larger transmittance values. In October, when ice is thinner and not much snow accumulated yet after summer melt, transmittance values are larger than 0.025. In situ measurements of light transmission before April-May are not available. However, in a recent compilation of under-ice light measurements, Katlein et al. (2019) report values of transmittance that do not exceed 0.01 in May. This is in agreement with the present estimations, except for the marginal areas. In these cases, though, the estimated transmittance values are still in the range of variability shown by Katlein et al. (2019) for the late springsummer period. Between late April and early June (i.e., before melt onset), reported values of light transmittance range between 0.0006 and 0.003 (see Supplementary Figure S8 in Oziel et al., 2019), in agreement with the present results of T < 0.01. Another newly published data set described in Castellani et al. (2020) and available on PANGAEA 1 shows a mode at transmittance values in May and June below 0.05, again in agreement with the estimates presented in this study.

Daily Evolutions From February to April 2012
As shown above, the use of the monthly means of each physical parameter gives a good evaluation of inter-annual changes in under-ice PAR on a monthly basis. However, for the timing of a bloom, daily estimates are preferred. Here we evaluate the influence of daily variations in snow depth and albedo on the under-ice PAR, but fixing the SIT to the monthly mean values and following the daily evolution of snow depth, incoming solar radiation and surface albedo. In particular, we are interested in when the daily PAR value increases above the threshold range of 0.5-5 µmol m −2 s −1 . In other words, when plankton is 1 https://doi.pangaea.de/10.1594/PANGAEA.902056 susceptible to be more active. Such a light threshold value is purely indicative, since a light threshold for net photosynthesis is not well observationally constrained or even well justified. Depending on the value of the threshold in the explored range and on the rate of change in irradiance, the alleviated onset of algal activity varies by days to a few weeks. Figure 10 shows the evolution of the under-ice PAR from 1 February to 30 April 2012, as well as the corresponding evolution of the snow depth and the albedo, for three randomly chosen locations in the Chukchi (72.97N/−171.29E), Beaufort (78.45N/−150.23E), and Laptev (77.09N/119.79E) seas. Note since the albedo and incoming solar radiation are centered around 14:00 GMT, these values represent the evolution at this time of day. Also shown is the algal growth range threshold (gray shading) based on Letelier et al. (2004) and Tremblay and Gagnon (2009). First, and as expected, we can observe the important dayto-day variations in amplitude of both snow depth and albedo and the influence these have on the under-ice PAR. Second, we notice the impact of the spatial variability of snow depth on the range of amplitude of the under-ice PAR between regions. As also shown in the monthly maps (Figure 9), less snow in the Laptev Sea compared to the Beaufort and Chukchi seas results in a substantial increase of the under-ice PAR earlier in the year. Finally, already by the end of February, the amount of under-ice PAR exceeds the lower threshold for initiating algal growth within the Chukchi Sea, though this is strongly controlled by the amount of snow on the ice and the albedo. In the Beaufort and Laptev sea locations, the lack of sunlight delays the timing of light under the ice until March, however, once there is sunlight, the under-ice PAR within the Laptev Sea location almost already exceeds the higher threshold for initiating algal bloom as the snow cover is considerably thinner than in the Beaufort Sea.

Impact of Homogeneous vs. Heterogeneous Snow Depth and Ice Thickness Distributions
Here we examine the impact of computing light transmittance through a heterogeneous snow depth and SIT distribution FIGURE 11 | Difference in monthly mean April PAR between using a heterogenous snow depth following Eq. 12 and assuming a mean ice thickness per grid cell versus a mean snow depth per satellite pixel. compared to assuming mean snow depth and ice thickness values per 25-km grid cell. For these comparisons, we use a mean SIT per grid cell to only examine the snow distribution impact, and vice versa, a mean snow depth to examine the influence of the SIT distribution. Further, only the month of April is shown as it is the month with the largest amount of light transmission through the ice prior to melt onset. While results are largely insensitive to year, years with thinner snowpack, or thinner ice will show overall larger effects of including a distribution versus mean values.
As expected, for both the ice thickness and snow depth distributions, the amount of under-ice PAR increases compared to using mean values per grid cell (Figures 11, 12). This is because thin ice and snow transmit much more light than thick ice and snow. With our parameterizations, the mean transmission is always higher than the transmission calculated with mean thickness -a result already obtained by Abraham et al. (2015) -in turn increasing the overall under-ice PAR for each grid cell. Nevertheless, the differences in under-ice PAR remain small, especially over the central Arctic where PAR increases just slightly, between 0 and ∼2 µmol m −2 s −1 (∼5% relative difference, or percentage of the average under-ice PAR it represents). Castellani et al. (2017) report a threshold value FIGURE 12 | Difference in April PAR between using a mean sea ice thickness per satellite pixel versus a fifteen ice thickness categories ITD heterogeneously distributed between 0 and 3h i , as described in section "Data Products," and assuming a mean snow depth per grid cell.
for algal bloom of 1.78 µmol photons m −2 s −1 which is at the upper boundary of the differences of PAR in the central Arctic. The effect of using mean values per grid cell thus does not have large effects on algal bloom, at least in the central Arctic. It is different for coastal regions in the Arctic basin and also in the Barents Sea, characterized by generally thinner snow cover as well as thinner ice, where differences in under ice PAR are large enough (i.e., >5 µmol m −2 s −1 , or >∼12.5% relative difference from the mean value) to impact the onset of algal bloom.
In regards to snow depth, increases in under-ice PAR above 4 µmol m −2 s −1 (10% of relative difference) are found mostly outside of the Arctic basin (i.e., Bering Sea, Baffin Bay, Barents Sea), but also at times within the Chukchi and Beaufort seas, as well as the Laptev and Kara seas that can reach 5 µmol m −2 s −1 (12.5% relative difference) or more. The largest increases in under-ice PAR within the Beaufort Sea are found in 2016, when there was relatively more transmittance of light as a result of less snow accumulation compared to other years (see also Figures 2, 4). Thus, the impact of using a snow depth distribution varies from year to year, depending on the average snow depth per pixel.
Similarly, for the ice thickness distribution versus mean SIT values we observe that larger differences on the order of at least 4 µmol m −2 s −1 occur when the sea ice is thinner, as for example close to the shelf but also in the Chukchi and Beaufort seas, especially during 2013, 2016, 2017, and 2018. As for Figure 10, April 2016 is the year when the differences between the SIT distribution and mean SIT per pixel are greater, especially in the Chukchi Sea with more than 5 µmol m −2 s −1 difference. This is due to the fact that for this year, the ice is overall thinner than for the other years. During some years and regions, the subgrid scale SIT distribution results in larger increases in under-ice PAR compared to sub-grid scale snow depth distribution (i.e., 2013 in Chukchi and Laptev seas), and in others the reverse is true (i.e., 2012 and 2015 in Beaufort Sea). Considering the threshold reported by Hancke et al. (2018), these sub-grid scale variations in snow depth and ice thickness may be important in terms of timing of under-ice algal growth and thus more research on how best to represent sub-grid scale snow depth and ice thickness distributions is warranted, especially during this transition from a multiyear to first-year dominated Arctic Ocean, with corresponding changes in surface roughness.

Impact of Snow Depth Trends on Transmission Through Snow
So far, our assessments have focused on the CryoSat-2 timeperiod and during the time of the year for which thickness observations are available. To assess longer term changes in light availability, not just in autumn and spring, but also in summer, we can examine how changes in snow depth and surface albedo are impacting the amount of light transmitted to the top of the ice (i.e., through the overlying snowpack). In particular, we are interested in the impacts of changing snow accumulation on light availability when light first becomes available and until the snow melts off the ice. Figure 13 shows trends from 1982 to 2018 in SnowModel-LG snow depths forced with MERRA-2 atmospheric reanalysis from February through July (see also Stroeve et al., 2020). Overall, winter and spring snow depth is found to be declining throughout the Arctic Ocean marginal seas, with slightly positive trends north of Greenland. In particular, statistically significant negative trends (at 95% confidence interval) are seen throughout the Beaufort, Chukchi, East Siberian, Laptev, Kara and Barents seas in February to May (between 2 to 8 cm per decade, depending on location), that increase to cover most of the Arctic Ocean also in June, and the central Arctic in July. Positive trends (∼1-3 cm per decade) north of Greenland are not statistically significant, but agree with positive trends in albedo in this region (Figure 14). Since no trends in precipitation are revealed in the reanalysis precipitation themselves (Barrett et al., 2019), the reduction in snow depth in the marginal seas is most likely a result of later freeze-up and earlier melt onset (e.g., Stroeve and Notz, 2018) that reduce the time over which snow can accumulate, rather than a change in precipitation. It also reflects a reduction in multiyear ice, as the longer an ice parcel can accumulate snow, the deeper the snowpack tends to be (e.g., Liston et al., 2020).
At the same time, and not surprisingly, the surface albedo also exhibits generally statistically significant negative trends over most of the Arctic basin in May, June, and July, reflecting earlier melt onset, earlier ice retreat (i.e., open water) and melt pond development. Negative albedo trends in February to April are mostly confined to the regions outside the Arctic basin (e.g., Baffin Bay, Barents and Bering seas). Negative albedo trends outside of the central Arctic prior to melt onset are largely the result of a lack of winter sea ice in more recent years, as the trends include both sea ice and open water regions. Small positive trends are also observed in February and March (∼0.02-0.06 per decade) over the areas with sunlight in the central Arctic, larger in February than in March. None of these trends are statistically significant, yet we do find they are consistent with slightly positive trends in snow accumulation that are statistically significant north of the Canadian Archipelago. On the other hand, positive albedo trends may also reflect problems in the albedo estimates during this time of year. As seen previously in Figure 10, albedo in 2012 is seen to reach nearly 1.0 in the Chukchi Sea in February which exceeds the albedo for new snow (0.90, Wiscombe and Warren, 1980). This likely points to a problem in properly accounting for the anisotropic reflectance of snow, especially under oblique solar and sensor zenith angles found at this time of year. Positive albedo trends (∼0.01 to 0.02 per decade) in April may also reflect changes in cloud cover. The APP-X data product produces an all sky albedo, whereby the amount of solar radiation underneath clouds is modeled based on estimated cloud top height, temperature and ice versus water clouds. Problems in accurately detecting cloud cover and cloud properties, combined with overall increases in springtime cloud cover (e.g., Wang et al., 2012) and cloud thickness (Huang et al., 2019) will influence the amount of incoming solar radiation as well as the surface albedo through modification of the spectral distribution of the incoming solar radiation (e.g., Grenfell and Maykut, 1977;Grenfell and Perovich, 2008).
The regions with negative snow depth trends, combined with negative albedo trends will play a significant role in increasing light availability under the ice. While we do not have ice thickness during that time, we can estimate how these snow depth changes, together with albedo changes, impact light transmission through the snow pack. These results are summarized in Figure 15. The largest positive trends in transmittance under the snow, exceeding 0.06 per decade are found in June in the Chukchi Sea and then everywhere there is sea ice in July. However, already in February small positive trends in under-snow transmittance are observed, that are statistically significant in the Chukchi Sea, in the East Greenland Sea and north of Novaya Zemyla (∼0.01 to 0.03 per decade). By April, statistically significant positive trends in under-snow transmittance occur throughout the marginal seas of the Arctic Ocean, reaching 0.04 per decade in the northern Barents Sea. This will cause a shift of the production season onset to earlier times in the spring in these regions. This will in turn affect the lower trophic levels, since there might be a mismatch between high food availability and events such as reproduction and spawning (Durant et al., 2007). Moreover, since the trend in snow cover is not uniform throughout the Arctic, this might lead to not only a temporal, but also a spatial shift of growth onset.

Uncertainty Analysis
As mentioned in section "Implementation, " we investigate the error of the under-ice PAR as a result of uncertainty on the FIGURE 13 | Snow depth trends (cm/yr) in February through July from 1982 to 2018 over Arctic sea ice as represented in SnowModel-LG  forced by MERRA-2 atmospheric reanalysis. Trends statistically significant at the 95% confidence level are indicated by black marks.
physical sea ice/snow parameters used to compute the light transmission as well as the values used in the radiative transfer ( Table 2). This evaluation was conducted for April of each year of the study and the resulting errors, averaged over the Arctic basin and for all 8 years (2011-2018) are summarized in Table 3. Overall, the uncertainty in the under-ice PAR from uncertainties in the physical sea ice or snow parameters ranges from ±1.5 to ±9.2 µmol m −2 s −1 , or a relative error of ∼4-22%. The uncertainty of ±20% on the SIT, gives the lowest errors in the under-ice PAR, while an uncertainty of ±10 cm in snow depth gives the highest PAR errors. This is followed closely by a ±5% uncertainty in SIC with a 7.2 µmol m −2 s −1 or (17%) and ±0.05 uncertainty in the surface albedo, which results in 5.0 µmol m −2 s −1 , or a relative error on the order of 12.2%. The larger impact of snow depth and albedo uncertainties relative to SIT on underice PAR agrees with Katlein et al. (2019), who found that ice thickness was an overall poor factor in determining overall light transmittance levels. Instead, snow cover was found to provide the primary limit on the amount of sunlight getting through the ice (i.e., less than 1% of light is transmitted below the ice in May), but light transmittance increases to 10% during the advanced phase of melt pond development, which is reflected in part by reductions in albedo during June and July (e.g., Figure 13). Uncertainties in ice concentration naturally allow for more (less) light absorption in open (closed) water areas.
On the other hand, uncertainties and relative errors for extreme values of the attenuation coefficients have lower overall relative errors compared to snow depth uncertainty, but errors can be as high as 10% for the large range in ice attenuation coefficients used. More important, however, are the large uncertainties in the under-ice PAR as a result of varying the i o,snow , contributing to an under-ice PAR uncertainty of 7.5 µmol m −2 s −1 , or relative error of 18.4%. In summary, uncertainties are within 20-30%, which is enough for a qualitative analysis of the pan-Arctic under ice light field, but not for a precise prediction of light levels in specific locations, nor for a precise prediction of timing for the bloom onset. Since the largest uncertainties are associated with snow properties (depth and i o ), this points to snow as the research area where most progress is needed.
Finally, the choice of SIC algorithm used can have a large impact in under-ice light transmission. While Kern et al. (2019) FIGURE 14 | Albedo trends per year in February through July from 1982 to 2018 over Arctic sea ice as derived from AVHRR. Note albedo trends are given for open water as well as ice-covered regions and thus large negative trends along the ice margins reflect in part earlier development of open water, while trends over the icepack in summer reflect also melt pond development. Trends statistically significant at the 95% confidence level are indicated by black marks.
found on a pan-Arctic scale differences between algorithms may cancel out, the NASA Team SIC underestimate those retrieved from visible imagery (e.g., from the NASA Moderate Resolution Imaging Spectroradiometer (MODIS) by 5-10%, whereas SICs from algorithms such as the Bootstrap algorithm (Comiso et al., 1997) overestimate MODIS-derived SICs by a similar amount. Since the Bootstrap algorithm has higher SICs than those from NASA Team, this will reduce the amount of light reaching the ocean surface, especially in the marginal ice zone where the differences between the two algorithms are most pronounced.

APPROPRIATENESS OF METHODOLOGY
In this contribution, we attempt to compute under-ice light levels based on current satellite products, and try to identify uncertainty sources. In this context, the two-level Beer-Lambert approach was the simplest available scheme, which takes into consideration the exponential decay of light through the snow/ice/water media, and gives reasonable transmittance retrievals, as compared with a significant number (N = 234) of data points from the seasonal ice zone (see Figure 1). The setting of AOPs, such as attenuation coefficients, were tuned using field observations over first-year ice, and may not necessarily apply over multiyear ice. However, we found that at least in regards to κ i , the value obtained of 1 m −1 differed little from classical multiyear ice derived values (1.5 m −1 ). These values are also in line with those of Grenfell and Maykut (1977) (see their Figure 5).
In terms of uncertainties in our under-ice irradiance calculations, the most important uncertainty source relates to non-existent or largely imprecise summer satellite products. In spring, when satellite products are less prone to error, uncertainties are smaller but still within 30% for each of the FIGURE 15 | Under-snow transmittance trends per year computed using the albedo and snow depth in February through July from 1982 to 2018 over Arctic sea ice. Trends statistically significant at the 95% confidence level are indicated by black marks. uncertainty sources. Next to surface albedo, snow in all its aspects appears as a dominant source of uncertainty. This includes the pixel-mean snow depth value, its sub-pixel distribution, and the description of radiative transfer within snow, which is underlined by the large errors on low transmittance values (see Figure 1), and also by the large impact on calculated under-ice PAR of the surface transmission parameter i o within snow (see Table 3). The new SIT product, derived from ice freeboards processed with the LARM algorithm and using snow data from SnowModel-LG improves the detection limit of CryoSat-2 for thin ice, with a lower freeboard detection limit of ∼2.5 cm (Landy et al., 2020). As evidence for this enhanced thin ice detection, we find that ice thicknesses < 0.5 m make up around 45% of measurements across the October SIT fields, and thicknesses <0.25 m make up around 31% of measurements. While the attenuation coefficients for ice based on first-year ice field data may not be the same as for multiyear ice, uncertainty in k i played a smaller role than the uncertainty in i o . These conclusions are in line with those of Katlein et al. (2019), which is based on a different observational dataset and with a slightly different approach.
To improve upon these uncertainties, calculations of underice light field would primarily benefit from more precise satellite products, spanning longer periods of the year (discussed further in "Looking to the Future"), in particular regarding ice thickness, snow depth and surface albedo. The sub-pixel distributions of ice characteristics are beyond reach for current observing systems. Yet these can at least be parameterized. For instance, Abraham et al. (2015) provide analytical calculation techniques to account for the sub-pixel snow depth distribution in underice light calculations. These can readily be used and at low cost. However, more research is needed to better understand what ice thickness distribution should be used as the ice cover has thinned and transitioned to a first-year ice regime. Generation of open water/leads, melt ponds and thick ice from ridging may differ from earlier studies, and in different regions of the Arctic, which will not only influence the ice thickness distribution, but also snow accumulation. Thus it remains unclear if a pan-Arctic representation of ice thickness and snow depth distribution is applicable as used in this study.
Regarding the formulation of the radiative transfer approach used, an easy recommendation would be to better constrain optical snow parameters from observations over various ice types, based on more observations, more precise, including more snow parameters. Another recommendation would be to use better representations of radiative transfer within snow and ice. Indeed, an easy objection to our methodology is that Beer-Lambert Law is overly simple and does not fully hold in snow, a highly scattering medium (Perovich et al., 2017). Another critique is that the surface scattering layer could be not as well defined in snow as it is in sea ice. While many studies have used the Beer-Lambert Law assuming an exponential attenuation of radiation with length of the medium (e.g., Maykut and Untersteiner, 1971;Bitz and Lipscomb, 1999;Castellani et al., 2017;Vancoppenolle and Tedesco, 2017), more physically elaborated approaches for radiative transfer in snow and ice do exist, in particular twostream methods, such as such as the delta-Eddington scheme of Briegleb and Light (2007), which relies on the specification of intrinsic optical properties of snow and ice, namely scattering and absorption coefficients. Two-stream schemes have better theoretical foundations and enable the representation of a wider range of optical phenomena (Dang et al., 2019). Yet such schemes are also more complex, more expensive and require information on the vertical profiles of temperature, salinity, density and microstructure (gas, brine, minerals, and grain size), information we do not have from satellites. Hence, at this stage, the use of more elaborated approaches does not warrant lower uncertainties in under-ice light distribution.
Finally, the impact of impurities on transmitted light intensity could also be considered. In particular, ice algae may reduce transmitted light when sufficiently abundant. Here we decided not to take them into account, because we found their effect to be generally low as compared with uncertainties on under-ice transmitted light, and largely uncertain because of unresolved spatio-temporal variations in ice algal content. Significant effects may occur where and when algae are ∼30 mg chl m −2 , and when snow is thin, which may occur in the Bering Sea, or during the season of algal bloom (generally after the month of April). One prerequisite before accounting for ice algae, would be to better constrain their abundance and spatial variability on scales relevant for our study. Massicotte et al. (2019b) show that to properly estimate primary production in the Arctic, we need to assess the spatial variability of sea-ice properties. By under-ice profiling platform (e.g., ROVs and towed nets) we can assess variability on a floe scale, but it is only by means of satellite products that we can provide important information on a pan-Arctic context. This study allows the characterization of the light field over an area large enough to capture the irradiance variability on a pan-Arctic scale. Furthermore, for the first time we are able to provide pan-Arctic maps of under-ice light at a spatial coverage and resolution of most sea-ice numerical models. Most large-scale sea-ice -ocean circulation models adopt the two layers Beer-Lambert approach to describe the radiative transfer, and they employ a SIT and snow depth distribution as in section "Impact of Homogeneous vs. Heterogeneous Snow Depth and Ice Thickness Distributions" to account for sub-grid variability. The present study is a first step in providing products that can be used to directly compare model outputs on a pan-Arctic scale. Moreover, for the first time, we quantify the effects of the sub-grid distribution on the under-ice light levels.

LOOKING TO THE FUTURE
Ideally, we would want to additionally map the underice light field from February through October on a pan-Arctic scale, accounting for SIT changes in summer, advanced snow metamorphism, development of melt ponds, and properly account for snow accumulation around ridges. Unfortunately, we lack summer SIT data and melt pond information to extend the characterization of under-ice irradiance beyond April. Temporally, we have been limited by the availability of sea ice freeboard and thickness products. CryoSat-2 provides profiles of SIT observations across the Arctic on a daily basis, but only at the satellite sub-cycle of 30 days does it provide complete pan-Arctic coverage and current algorithms using radar altimetry remain limited to October through April. Since the change from complete darkness to enough sunlight for photosynthesis to occur can happen over a few days, daily coverage is needed. Blending of CryoSat-2 satellite retrievals with Sentinel 3A and 3B, can increase temporal resolution to 10-day frequency (Lawrence et al., 2019). Nevertheless, prior to the melt season, snow depth is arguably more important in limiting light availability and timing of ice algal blooms.
In summer, the most important data gaps to fill are the ice thickness and melt ponds. ICESat-2 offers the ability to provide sea ice freeboards in summer, yet the conversion to total SIT is difficult as this requires information on ice density which is highly variable during the summer melt season. Instead, one could consider using the ICESat-2 freeboards to represent the height of the strongly scattering surface layer (SSL). Typically, the layer has geometric depth 3 -5 cm and effective scattering about two orders of magnitude larger than the interior ice. The ice between the SSL and sea level ('drained layer' DL) also has enhanced light scattering properties, typically one order of magnitude above that of interior ice (see Light et al., 2015, Table 1). Despite their small vertical extent, but because of their strong scattering, these two layers are estimated to have optical depth approximately two orders of magnitude larger than the entire portion of the ice that sits below sea level. Because light penetration will depend most strongly on the portions of the ice with the largest optical depth, it may be possible to estimate transmittance based on the thicknesses of the SSL and DL, as derived directly from summer freeboards. The assumption is that the optical depth of ice above freeboard is about many times more than the optical depth of ice below freeboard, and therefore the amount of ice floating above the water is more important for determining how much light reaches the ocean beneath the ice than the actual SIT.
Once the melt season starts and the snow disappears, melt ponds dominate the amount of light entering the upper ocean beneath the ice. Melt ponds develop in early summer of sizes typically < 10 m, and evolve into a temporally and spatially heterogeneous network of melt ponds intermixed by bare ice and snow-covered ice regions. While melt ponds are not difficult to detect using high resolution satellite imagery, such as that from meter-scale WorldView optical imagery (e.g., Wright and Polashenski, 2018), pan-Arctic coverage and long time-series are not available from such high resolution imagery. Instead, moderate resolution satellite sensors, such as the NASA Moderate Resolution Imaging Spectroradiometer (MODIS) satellite have been used to map melt ponds at a pan-Arctic scale (e.g., Rösel and Kaleschke, 2012). MODIS has a spatial resolution of 250 m for the first two visible channels and 500 m for the other five visible channels. Since the melt ponds have a relatively small size in relationship to the large footprint of the satellite image, several ice types are present within each satellite pixel, requiring some sort of spectral unmixing to retrieve the melt pond fraction. Spectral mixture analysis (SMA) (e.g., Tschudi et al., 2008;Rösel and Kaleschke, 2012) and Multiple Endmember Spectral Mixture Analysis (MESMA) (e.g., Yackel et al., 2018) have been used with limited success to map sub-pixel fractional area of surface types, yet reference spectra for ice types can vary considerably and cause large errors in melt pond fraction retrievals. Currently there are no up-to-date pan-Arctic melt pond products being provided by the science community.
Similarly, as the ice cover has transitioned from one dominated by multiyear ice year-round to one dominated by firstyear ice, the dynamics of the ice cover have increased, creating more occurrence of ice fracture and lead formation (e.g., Rampal et al., 2009). Leads result in new ice formation, and thus thin ice with little to no snow cover, and thus play an important role in light transmission before the melt season starts. Evidence from an Arctic fjord found that 77-86% of incidence light from 0.4 to 0.7 µm was transmitted through leads (Taskjelle et al., 2016). Changes in albedo may in part reflect increased fracturing, leads and thin ice formation. On the other hand, evaluation of the APP-X data set shows non-physical artifacts in albedo that have an impact on the under-ice PAR calculations, and more improvements are needed to develop a long-term and reliable surface albedo product.
Finally, ridges play an important role in redistributing snow accumulation. While several approaches exist to quantify surface roughness from satellite, including laser altimetry (e.g., Landy et al., 2015;Petty et al., 2016), multi-angular optical imagery (e.g., Nolin and Mar, 2019), and SAR backscatter (e.g., Fors et al., 2016), the small footprint of the ICESat-2 ATLAS photoncounting instrument (∼15 m), together with dense along-track spacing (∼70 cm) and precise elevation measurements (<10 cm) will allow for mapping of ridge heights. This is because ICESat-2 will provide consecutive elevation measurements to define the ridge height and edges with sufficient precision. This information could be used in conjunction with SnowModel to map snow depth at high spatial resolution and then downscale to the pan-Arctic scale.

CONCLUSION
Arctic marine ecosystems have adapted to extremes in light conditions, taking advantage of the short time-period over which primary production can occur. Since satellites sensors cannot 'see through' the sea ice, the response of Arctic primary production to changes in sea ice and snow properties can only be determined with continuous in situ monitoring of algae and phytoplankton stocks, together with the main physical drivers affecting their phenology: under-ice light field, sea ice and snow properties, temperature and salinity of the ocean surface, nutrients and mixed layer depth. However, the lack of long-term data covering different regions and seasons makes this challenging at present. In this paper we have attempted to map one component important for ocean primary production, the light under the ice on a pan-Arctic scale using currently available satellite-based products on snow and ice conditions together with parameterizations established from cruises and other in situ observations. Our estimates of under-ice PAR are based on the Beer-Lambert Law, and our interactions of light with snow and ice are based on empirical parameterizations from in situ data collections that depend at the moment solely on the ice thickness and snow depth, whether or not the snowpack is melting, as well as the surface albedo. Key uncertainties in our under-ice PAR calculations relate to how precise current parameterizations are, how best to distribute the mean snow and ice thickness within a 25-km grid scale, as well as errors in satellite-retrieved variables of surface albedo, ice thickness and concentration, and how best to model snow accumulation.
During springtime, the high albedo of snow and its capacity to attenuate solar radiation means that snow depth variability largely controls under-ice PAR variability from year-to year, with some years showing increased light penetration especially in the Beaufort Sea in April. During October, ice thickness, albedo and lead fraction play a more important role as there is little snow on the ice. Ice algal in particular have adapted to very low light levels (e.g., Berge et al., 2015), with the first bloom happening in the bottom layer of the sea ice already in spring when light for photosynthesis becomes available. Thus, despite little light entering the Arctic Ocean in February and March, inter-annual snow depth variability can lead to more light availability even during these months that is enough to initiate algal growth, whereas long-term declines in snow depth over time may be shifting the timing of under ice algal blooms to earlier in the year. This has important implications on marine ecosystems, as ice algae are an early food source for certain pelagic grazers (e.g., Leu et al., 2015).
New satellite sensors, such as ICESat-2, could enable these under-ice maps to be extended to be year round, whilst ongoing in situ field programs, such as EcoLight and MOSAiC, will allow us to refine the algorithms and more accurately predict light penetration through a range of snow and ice conditions. While there are significant challenges to overcome before we routinely produce daily pan-Arctic under-ice light maps using satellite observations, we feel new satellite-derived products and focused in situ campaigns are closing the gap.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: ftp://sidads.colorado. edu/pub/projects/SIPN/Ecolight. All AVHRR Polar Pathfinder albedo and incoming solar radiation data are available from NOAA National Centers for Environmental Information (https: //www.ncdc.noaa.gov/cdr/atmospheric/extended-avhrr-polarpathfinder-app-x). Sea ice concentration data is available from the National Snow and Ice Data Center Sea Ice Index (https: //nsidc.org/data/seaice_index).