Seasonal and interannual variations in the propagation of photosynthetically available radiation through the Arctic atmosphere

The Arctic atmosphere–surface system transmits visible light from the Sun to the ocean, determining the annual cycle of light available to microalgae.This light is referred to as photosynthetically available radiation (PAR). A known consequence of Arctic warming is the change at the atmosphere–ocean interface (longer ice-free season, younger ice), implying an increase in the percentage of PAR being transferred to the water. However, much less is known about the recent changes in how much PAR is being transferred by the overlaying atmosphere. We studied the transfer of PAR through the atmosphere between May 21 and July 23 at a pan-Arctic scale for the period ranging from 2 000 to 2 0 16. By combining a large data set of atmospheric and surface conditions into a radiative transfer model, we computed the percentage of PAR transferred to the surface. We found that typical Arctic atmospheres convey between 6 0 % and 7 0 % of the incident PAR received from the Sun, meaning the Arctic atmosphere typically transmits more light than most sea ice surfaces, with the exception of mature melt ponds. We also found that the transfer of PAR through the atmosphere decreased at a rate of 2.3% per decade over the studied period, due to the increase in cloudiness and the weaker radiative interaction between the atmosphere and the surface. Further investigation is required to address how, in the warmer Arctic climate, this negative trend would compensate for the increased surface transmittance and its consequences on marine productivity.


Introduction
Earth's obliquity gives rise to large seasonal variations in the input of solar irradiance to Arctic marine ecosystems. Notably, it limits primary productivity during several months of the year, particularly during the polar night. During other times of the year, sea ice (including overlying snow) is what most constrains the amount of solar irradiance used by algae, named photosynthetically available radiation (PAR), in the water column. The impacts of ongoing climate-related changes on sea ice transparency to PAR have been emphasized in several recent studies (e.g., Arrigo et al., 2012;Assmy et al., 2017;Horvat et al., 2017). Significant reduction in the transfer of PAR through the atmosphere has also been observed at the pan-Arctic scale, impacting the marine primary production . PAR estimations at the sea surface by Bélanger et al. (2013), however, ignored the radiative interaction (i.e., multiple back and forth diffuse reflections) between the atmosphere and the highly reflective surface in the presence of sea ice. Laliberté et al. (2016) included this interaction in their PAR model and, based on in situ observations, showed its significance. Here, we further examined the variations in the propagation of PAR through the atmosphere over the Arctic Ocean while taking into account the interactions between the atmosphere and surface. We focus on the spring-to-summer transition when the solar irradiance is highest and both the ocean surface and the atmosphere change rapidly due to rising air temperature and melting snow and sea ice.
In this study, the surface is a horizontal layer that includes only the air-ocean interface for open waters and that extends enough vertically to encompass sea ice when it is present (including the air-ice and ice-ocean interfaces). Following this definition, we define the surface transmittance for PAR, T s (%), as the ratio of downward irradiance integrated over the PAR spectral range (400-700 nm) in the water just below the sea surface or the bottom of the ice, E d (PAR, 0 -) (moles photons m -2 d -1 ), and downward irradiance in the air just above the sea surface or the top of the ice, E d (PAR, 0 þ ) (moles photons m -2 d -1 ): Accordingly, the percentage of PAR transferred through the atmosphere, T a (%), is defined as the ratio between E d (PAR, 0 þ ) and downward irradiance at the top-of-theatmosphere, E d (PAR, TOA) (moles photons m -2 d -1 ): Equation 2 presents the downward irradiance at the surface normalized by its top-of-the-atmosphere value. Note that T a , just like T s , could be seen as the "atmospheric transmittance," except that it is significantly impacted by the boundary conditions. In other words, the radiative properties of the bottom layer at the base of the atmosphere are much more variable than those for T s and can influence the amount of light available at the surface. For this reason, we will refer to the "percentage of PAR transferred through the atmosphere" or "normalized downwelling irradiance" rather than "atmospheric transmittance" to avoid confusion with the latter concept. Finally, in the chosen framework, E d (PAR, 0 -) is available to Arctic marine primary producers and results from E d (PAR, TOA) and the balance of T s and T a : T s variations in the Arctic Ocean are mostly driven by changes in sea ice concentration and physical properties, including thickness, brine and air content, presence of melt water, and the state of the snow cover (Perovich, 2005;Arndt and Nicolaus, 2014). Because of their effect on E d (PAR, 0 þ ), the most critical factors controlling variations in T a are cloudiness (Rabbette and Pilewskie, 2002;Bélanger et al., 2013) and surface albedo (Grenfell and Perovich, 2004).
Sea ice cover and thickness are decreasing (Haas et al., 2008;Stroeve et al., 2012), while the older ice is replaced by more dynamic (more leads) and younger ice (larger melt pond coverage; Comiso, 2012;Nicolaus et al., 2012;Wang et al., 2016), resulting in an increase of T s . In parallel, a general increase in cloudiness (Eastman and Warren, 2010) is especially noticeable in the spring and summer (Wang and Key, 2003), resulting in a decrease in T a . This work focuses on evaluating the transfer of visible radiation in the atmosphere, T a , and its change in space and time over the Arctic Ocean. We decomposed the Arctic Ocean into three regions, defined by their ice cover, and the marine productive season into four periods of the year, relevant to the spring bloom and centered on the summer solstice. Using satellite remote sensing data from multiple sources combined into a radiative transfer model, we computed T a over 17 years for each region and period of the year. Interannual and seasonal trends are calculated and interpreted. Finally, we verify how T a compares in magnitude to values of T s .

Model
E d (PAR, 0 þ ) and E d (PAR, TOA) were computed every day between May 21 and July 23 for every 1 Â 1 pixel above the Arctic circle (66 N) for nonland areas between 2000 and 2016, resulting in 8,640 pixels Â 64 days Â 2 estimates per year.
To compute the propagation of PAR through the Arctic atmosphere, we used the Santa Barbara Discrete ordinate Atmospheric Radiative Transfer (SBDART) computer program (Ricchiazzi et al., 1998), a one-dimensional atmospheric radiative transfer model. The input variables were cloud fraction, cloud optical thickness, cloud top height, water vapor amount, total ozone amount, aerosol optical depth (AOD), and the surface albedo (see below). SBDART has been used extensively to propagate irradiance through the atmosphere (e.g., Shupe and Intrieri, 2004;Su et al., 2007;Grenfell and Perovich, 2008;Painter et al., 2012). For each pixel, two SBDART runs were made every day, one for the clear part of the sky, the other for the cloudy part. They were then blended conservatively according to the cloud fraction. Each run yielded spectrally resolved downward irradiance at a 3-nm spectral resolution and expressed in W m -2 mm -1 . The instantaneous estimates were then converted into daily values, assuming constant atmospheric and surface conditions throughout the day for the temporal integration (Frouin et al., 2003). Spectral results were integrated over the 400-700 nm band to achieve PAR estimates.

Inputs
All inputs from satellite data discussed in this section are summarized in Table 1, where the products, associated references, spatial/temporal resolutions, time span, and website of the satellite-derived inputs used to compute the above surface downward planar irradiance are presented.
An atmospheric data set was built from the Moderate Resolution Imaging Spectroradiometer (MODIS) Level 3 Atmosphere Gridded Product (MOD08 and MYD08; Platnick et al., 2017), which contains the diurnal (daytime) cloud fraction, cloud optical thickness, cloud top height, water vapor and ozone content, and AOD (10.5067/MODIS/ MYD06_L2.006, 10.5067/MODIS/MOD06_L2.006). These atmospheric properties were retrieved from multiple observations collected during different times of the day and averaged on an equal-angle global grid (Hubanks et al., 2019).We downloaded data from both Terra and Aqua platforms and averaged them to obtain daily fields on a pixel basis.
We generated a surface albedo following a large-scale generic model of the seasonal evolution of albedo for the visible spectral range, modified from Perovich and Polashenski (2012) and Perovich et al. (2011;see also Perovich et al., 2002;Perovich et al., 2007;Light et al., 2015a). The visible sea ice albedo values were approximated from the spectral albedo measurements of Perovich (1996) and Light et al. (2008Light et al. ( , 2015b). This empirical model distinguishes six transitions in the state of sea ice during a year ( Figure 1). It takes into account the sea ice age, sea ice concentration, and a number of melting or freezing phases that are detected from space using passive microwaves. The phases were determined from the dates corresponding to early melt onset, melt onset, early freeze onset, and freeze onset (Markus et al., 2009). These onsets are based on surface dielectric properties changes, Art. 9(1) page 2 of 17 Laliberté et al: Changes in the propagation of PAR through the Arctic atmosphere monitored using the microwave signature of the surface. After winter, the early melt onset is the date when the sea ice surface layer emissivity first increases sharply due to the change in volumetric moisture content. The surface layer then undergoes melt-freeze cycles for some time, but eventually reaches the melt onset, a state when sea ice persistently stays under melt conditions, with signature approaching that of a blackbody in the microwave domain. The two other markers, the early freeze onset and freeze onset go in the opposite direction. The early freeze onset date is determined by the microwave signature showing evidence that the surface is starting to refreeze, and then freeze onset date is determined when the surface emissivity expresses persistent cold winter conditions.
The surface albedo was determined assuming specific sea ice albedo for each phase. When a pixel was icecovered, its albedo progressed following a first-year ice (sea ice age indicating that it has not survived a melt period) or perennial ice scenario (sea ice age indicating it has survived a melt period; Tschudi et al., 2016). Although determining whether a parcel of first-year ice that has not survived a melt period was going to disappear over the summer was not possible, this unknown had little impact once albedo was combined with the sea ice concentration and melt/freeze indicators (Perovich et al., 2011). For both scenarios, a high stable albedo represented the cold winter phase. The first transition was triggered when the early melt onset date occured, leading to  a small albedo decrease in the PAR spectral range. Then, after the melt onset was detected, a slow linear decrease was applied over 15 days to reproduce the snowmelt period. After that, a faster albedo linear decrease was used over a shorter period of 6 days to represent the melt pond formation. Subsequently, we represented the pond evolution by a small and constant decrease in the albedo, lasting until the arrival of the early freeze onset. The albedo increased slowly after the early freeze onset and faster after the freeze onset date, leveling at its initial winter value. At last, the surface albedo (A s ) varied as a function of time and was a weighted mean of the fraction of sea ice and open waters (Cavalieri et al., 1996) with a variable sea ice albedo (A i ) and a constant water albedo (A w ) of 0.08 (roughly representative of sun at 60 zenith angle under a clear sky). Accordingly, we have where the sea ice albedo is allowed to vary with time. When the sea ice concentration was available, but the early melt onset, melt onset, early freeze onset, or freeze onset were not, A s was computed using an A i of 0.40. This situation happened mainly in regions with low ice concentration, where we expect the existing ice to be free of snow. We did not compute E d (PAR, 0 þ ) when cloud fraction, cloud optical thickness, or surface albedo were not available, which generally occurred over land-dominated pixels or near the North Pole where satellite orbit inclination limits data availability. When cloud fraction, cloud optical thickness and surface albedo were available, but cloud top height, ozone, water vapor content or AOD were not, missing values of the latter were replaced by climatological values, as these geophysical variables are not the main drivers of variations in E d (PAR, 0 þ ) in the Arctic (Barrientos Velasco et al., 2020). We calculated the daily climatology on a pixel basis from all available values over the 17-year MODIS data set. In the case of aerosol optical thickness, data were usually only available in summer because the presence of clouds and sea ice at other times prevented its retrieval (Tomasi et al., 2015). If satellite retrievals and even climatological values were nonexistent, we filled the gaps on the basis that AOD was 0.15 at the beginning of May and assumed to decrease linearly to 0.09 at the end of July (Tomasi et al., 2007;Breider et al., 2014). Finally, passive sensors do not allow to characterize the vertical distribution of clouds, but we found that changing the altitude of a given cloud deck in radiative transfer simulations generally resulted in E d (PAR, 0 þ ) changes of small magnitude (below 3%).

Comparison with a standard algorithm
Comparing the resulting E d (PAR, 0 þ ) product with near contemporaneous in situ measurement would not provide much insight because of the unknown subpixel spatial heterogeneity and the large spatial extent of the product. Comparing our product with other existing reanalysis products would also be difficult because of the lack of operational E d (PAR, 0 þ ) products. Besides, readily available surface radiative flux products (e.g., Extended AVHRR Polar Pathfinder data set; Key et al., 2016) are broadband, and their conversion to a PAR range is highly dependent on atmospheric conditions and solar elevation (e.g., Frouin and Pinker, 1995). However, an operational E d (PAR, 0 þ ) algorithm used for primary production assessment exists and was compared with our product for its evaluation. We compared the E d (PAR, 0 þ ) obtained with our method (hereafter referred to as SBDART) to the standard E d (PAR, 0 þ ) product of the NASA Ocean Biology Processing Group (OBPG; Frouin et al., 2003), which has been successfully evaluated in the Arctic (Laliberté et al., 2016). The OBPG algorithm computes E d (PAR, 0 þ ) based on an energy budget approach, assuming that clouds are decoupled from the atmosphere and that the photons that are neither absorbed nor scattered out of the atmosphere reach the ocean. The OBPG E d (PAR, 0 þ ) product is currently limited to open waters because this approach does not address (1) the losses of photons due to absorption in sea ice, (2) the strong coupling between the surface and atmosphere over sea ice, and (3) the surface reflectance anisotropy of the changing sea ice surface (see Laliberté et al., 2016). Here, the OBPG E d (PAR, 0 þ ) product was downloaded for both Terra and Aqua and subsequently averaged. All data from the two methods were matched temporally and compared on a pixel-by-pixel basis. We computed linear Type II regression major axis using the R package lmodel2 (Legendre, 2018) along with Pearson's correlation coefficient and the root-mean-square difference (RMSD):

Analysis
To assess the changes in Arctic T a when substantial marine primary production is expected (Ardyna et al., 2013), we divided the Arctic Ocean into three regions and four periods. The regions were defined as a function of their ice regime above 66 of latitude. This approach is justified by the coupling between the surface type cover and the atmospheric properties (Yu et al., 2019). Following Yu et al. (2019), the permanent ocean region was associated with all grid elements for which the sea ice concentration was consistently below 15% over 99% of all four periods considered for this study. In contrast, the permanent ice region was associated with all grid elements for which the sea ice concentration was above 15% over 99% of the time. The transient ice region was set to all grid elements not falling in the previous two categories. The permanent ocean region consisted mostly of the Greenland, Norwegian, and Barents seas. The transient ice included the Kara, Laptev, and East Siberian seas, as well as the Chukchi and Beaufort seas and Baffin Bay. The permanent ice region covered the Central Arctic Ocean and extends toward the West, as well as through the Canadian Archipelago. This division of the three regions is presented in Figure 2. Temporally, the productive period was divided into spring, late spring, early summer, and summer, each lasting 16 days (May 21 to June 5, June 6 to June 21, June 22 to July Art. 9(1) page 4 of 17 Laliberté et al: Changes in the propagation of PAR through the Arctic atmosphere 7, and July 8 to July 23). By doing so, we covered the timing when the maximum ice-covered and ice-free Arctic primary production is expected. We first examined the seasonal evolution of E d (PAR, 0 þ ) and the time of year when it reached its maximum value. Because different pixels covered different effective sizes of the surface, we used the area-weighted operator to bin the variables. The area-weighted median of E d (PAR, 0 þ ) seasonal time series were examined separately for each region and year. On the same time series, we also identified the day of annual maximum by applying a 7-day moving average. A longer window would have blurred the variability, and a shorter window was susceptible to establishing the maximum to short-lived events such as isolated clear-sky days. Because the evolution of E d (PAR, 0 þ ) varies greatly due to the seasonal insolation cycle, we also looked at its normalized value (T a ; Equation 2), allowing greater focus on the role of atmosphere and surface properties. Two aspects of T a were studied: the seasonal cycle and the multiyear trend. We derived the seasonal evolution for each region by calculating the median T a of areaweighted pixel value on each day of the year. The seasonal evolution was derived for each individual region over all years. We used the median to reduce the gap created by missing MODIS-Terra data in the second half of June 2001 (MODIS-Aqua was only launched the following year). This analysis allowed us to obtain an overview of the seasonal evolution of T a . We derived the multiyear trend for each region and for each period of the year (6 Â 4 ¼ 24 trends). The area-weighted median of each region period and each year was computed, and a linear least-squares regression was fitted on the 24 multiyear trends from 2000 to 2016. This approach allowed us to verify if T a had changed over recent times.

Results and discussion
3.1. Inputs Figure 3 provides a general description of the seven geophysical variables used as input in the SBDART computer program to determine the percentage of PAR transferred through the atmosphere, T a , over each of the 8,640 grid cells for 64 days each year between 2000 and 2016. The median, interquartile range, and 10th and 90th percentiles of the distributions were computed from the satellitederived variables and the albedo model (for numerical values, see Table S1), and the retrievals contributing to a single MODIS daily average were binned from multiple overlapping orbits (Hubanks et al., 2019). Consequently, the values presented in Table S1 represent a large amount of data (and multiple environmental conditions) and can be considered as climatological values representative of the Arctic. Cloud fraction was generally high over the Arctic around the summer solstice, with a moderate cloud optical thickness and low cloud altitudes, reflecting the prevalence of low-level clouds (Shupe and Intrieri, 2004). There was a higher cloud fraction over open waters. Visible surface albedo was highly heterogeneous through space and time for the regions with sea ice and decreased with time of year. The Arctic atmosphere was generally dry, especially in spring, but became more humid as firstyear sea ice melted and leads developed (Curry et al., 1996;Kay and Gettelman, 2009), or as temperature  Arctic AOD median was 0.13, which is somewhat higher than when only valid retrievals from MODIS are considered (median of 0.10). The homogeneity of AOD values reported in Table S1 is due to the fact that 80% of the estimates were from the climatology. Fortunately, the above surface PAR is not sensitive to minor variations in AOD. For example, in Baffin Bay (74 , -74 ) on June 1 under a clear sky or cloudy sky and wet snow on top of sea ice, reducing an AOD of 0.1 by 30% led to a reduction of above surface PAR of <1%.

Comparison with a standard algorithm
We compared the E d (PAR, 0 þ ) estimates obtained with the SBDART and OBPG methods. The latter was limited to open waters, defined as pixels having sea ice concentration lower than 10%. A total of 118,996 matching pixels yielded a regression slope of 0.72 and an intercept of 13.88 with a Pearson's correlation coefficient of 0.86 and an RMSD of 7.62 moles photons m -2 d -1 (Figure 4a). In general, OBPG PAR estimates were higher than the SBDART estimates over open waters. The two methods employ largely different approaches, making a step-bystep comparison difficult. However, we speculate that our prescribed constant and conservative open water albedo value versus the estimate of the instantaneous surface albedo value in the OBPG method may cause a bias between methods. Frouin et al. (2018) used the aerosol and cloud properties from MERRA-2 hourly data to assess the OBPG algorithm uncertainties and also found that the OBPG method was slightly overestimating PAR. A weak overestimation by the OBPG algorithm was also observed when compared to in situ data collected shipboard from the Arctic (Laliberté et al., 2016) or from a moored buoy on the western side of Canada (Frouin et al., 2003). Nevertheless, the models are in good agreement, with a relatively low RMSD and a high correlation coefficient. Spatial patterns in OBPG and SBDART products ( Figure 4b) are consistent, but a definitive advantage of using the SBDART method is that PAR can also be estimated in all sea ice conditions.

Seasonal evolution of E d (PAR, 0 þ )
The Arctic undergoes a strong seasonal cycle of E d (PAR, 0 þ ) due to the seasonal insolation cycle, with at least 1 day (summer solstice) when the sun is above the horizon for 24 h in the study period. Accordingly, the day of the summer solstice could be expected to be the day with the greatest incident irradiance, but in fact that is rarely the case. In this section, we examine the seasonal variations in E d (PAR, TOA), E d (PAR, 0 þ ) under clear skies and actual sky conditions, and the distribution of dates for which E d (PAR, 0 þ ) reached its maximum every year ( Figure 5).
As a reference, the E d (PAR, TOA) peaked on the day of summer solstice and exhibited symmetrical decrease at decreasing and increasing days of the year ( Figure 5). The E d (PAR, 0 þ ) calculated for a clear sky scenario also peaked on the day of summer solstice and showed a symmetrical distribution with time for the permanent ocean region where surface albedo was distinctively stable over all periods (Table S1). For the two other regions where the surface albedo experienced a considerable decrease between May 21 and July 23, the E d (PAR, 0 þ ) variation pattern under clear sky was slightly asymmetrical and peaked earlier (Figure 5), which suggests coupling between the atmosphere and the surface . The E d (PAR, 0 þ ) calculated when accounting for actual sky conditions was strongly asymmetrical and peaked up to 4 weeks before the summer solstice, with an average of 10.5 days before the summer solstice. This early peaking is consistent with what was reported by Bernhard et al. (2007), who observed a day of annual maximum occurring 3 weeks before the summer solstice using their ground-based spectroradiometer located 300 m from the Chukchi Sea in Utqiagvik (Barrow), Alaska. Over the studied years, the transient ice region had the earliest day of annual maximum, with a median of 16 days before the summer solstice. The region with the latest annual E d (PAR, 0 þ ) maximum was the permanent ocean region with a median value of 7 days before the summer solstice. The intuitive notion that the date of the summer solstice corresponds to an annual maximum in E d (PAR, 0 þ ) is not accurate in the Arctic. In other words, the seasonal insolation cycle is indeed a central factor determining the maximum observed E d (PAR, 0 þ ), but the annual cycle of Arctic environmental factors as well as synoptic events also have a significant effect. In fact, regions with similar top-of-the-atmosphere incident solar irradiances (i.e., in Figure 2, permanent ocean and transient ice regions cover similar latitudes) exhibit very different E d (PAR, 0 þ ) values and great interannual variability ( Figure 5). For the rest of this article, we normalized E d (PAR, 0 þ ) by E d (PAR, TOA) to obtain the percentage of PAR transferred through the atmosphere (T a ; Equation 2, also named normalized downwelling irradiance) and inquire further into this variability.

Seasonal evolution of T a
Variations in T a are expected to be much less dependent on the insolation than E d (PAR, 0 þ ). Nonetheless, there was a small, second-order, remaining dependency due to the more prevalent slant path of the incoming sunlight through the atmosphere away from the summer solstice. As we restricted our study to a brief period symmetrically positioned on either side of the summer solstice, we made sure that only an increase in cloudiness, a decrease in surface albedo, or a combination of both, can cause a transition of large magnitude in T a (Gardiner, 1987;Fitzpatrick et al., 2004;Taskjelle et al., 2017a). Figure  6 shows the median seasonal evolution (2000-2016) of T a found for every region. To explore the impact of clouds and sea ice on the propagation of PAR through the Arctic atmosphere, we computed T a for two additional scenarios. Keeping everything else constant, we computed the seasonal evolution for T a under clear skies (no clouds; T a clear sky in Figure 6) and for open water conditions (no sea ice; T a open ocean in Figure 6). In the figure, the lines appear rather smooth because of the spatial and temporal averaging, which blurs the day-to-day variability and emphasizes the bulk seasonal evolution of each region.
Relative to T a (red line in Figure 6, T a ), when removing the clouds from the simulations (turquoise line in Figure  6, T a clear sky), the percentage of PAR transferred through the atmosphere increased. The effect of clouds can be seen by the difference between the red line and the turquoise line in Figure 6. When sea ice is removed from the simulations (yellow line in Figure 6, T a open ocean), the normalized downwelling irradiance decreased (except on the right panel where the median T a is plotted over the median T a open ocean). The effect of sea ice is the difference between the red line and the yellow line. In all cases, the median percentage of PAR transferred through the atmosphere for the clear sky scenario remained above 85%, but always below 93%. This finding means that where the Arctic is cloud-free, we can expect that, in general, between 7% and 15% of E d (PAR, TOA) would be lost through the cloudless atmosphere before reaching the surface. When averaged over the four periods, T a clear sky was 91%, 90%, and 87% for the permanent ice, transient ice, and permanent ocean regions, respectively. When instead of removing clouds, we removed the sea ice (T a open ocean), the normalized downwelling irradiance was between 46% and 32%. When averaged over the four periods, we find T a open ocean of 39%, 42%, and 37% for the permanent ice, transient ice, and permanent ocean regions, respectively. These results allow us to further explore the general impact of clouds. Considering the difference between T a clear sky and T a open ocean, we may conclude that clouds are responsible on average for half (52%, 48%, and 50% for the permanent ice, transient ice, and permanent ocean regions, respectively) of the reduction of PAR being transferred through the atmosphere. Because, in reality, the Arctic surface is generally highly reflective, the normalized downwelling irradiance is higher than that. This fact is illustrated in both the regions of permanent ice and of transient ice, where early in the season (spring period), this reduction of normalized downwelling irradiance (due to clouds, gases, and aerosols) was effectively cut by one-half and one-third, respectively, due to the presence of reflective sea ice. The surface, coupled to the overlying atmosphere, strongly affects the visible radiation available at the surface. This analysis reveals the sensitivity of the normalized downwelling irradiance to the parametrization of albedo and underlines the importance of having a proper surface albedo parameterization in the context of Arctic primary production modeling, especially when considering leads and the marginal ice zone Blais et al., 2017).

Decadal trends in T a
For this section, we focus on the interannual changes in T a . Figure 7 presents the trends in T a from 2000 to 2016 separately for the three regions and four periods between May 21 and July 23. Overall, 10 of the 12 trends showed a decrease in T a with time, including five that were statistically significant (P < 0.05).
As seen on Figure 7, the Arctic atmosphere, in general, shifted from relatively transparent to more opaque to PAR over the years. In the spring (May 21 to June 5), we found negative slopes for all three Arctic regions, with a mean of -1.9%/decade and a statistically significant slope for the transient ice regime (P < 0.05). Late spring was also characterized by negative T a trends over all regions, with Laliberté et al: Changes in the propagation of PAR through the Arctic atmosphere Art. 9(1) page 9 of 17 a mean of -2.4%/decade and significant slopes for the transient ice. Early summer and summer show weaker trends with means of -0.9%/decade and -1.4%/decade, respectively. In summary, the permanent ice region displayed negative trends for the spring and late spring, but a positive trend in the early summer and no trend for the summer (none were significant). The transient ice region displayed negative and significant trends over all periods. In other words, the confidence intervals set at the conventional boundary suggest that there is a trend for the transient ice region over all periods. The permanent open ocean region showed negative trends for all periods, with the summer period as statistically significant, which, due to the definition of that region, points to an increase in cloudiness. Over the whole Arctic (considering all pixels weighted by their area), and for all considered periods (from May 21 to July 23), T a decreased on average by 2.3%/decade (statistically significant) between 2000 and 2016.
Since E d (PAR, TOA) is stable over time, trends in E d (PAR, 0 þ ) showed the same patterns as T a . The overall decrease in E d (PAR, 0 þ ) between 2000 and 2016 is consistent with the trends reported by Bélanger et al. (2013; their figure 1) for the 1998-2009 period. Interestingly, when computing the equivalent of Figure 7 for the open ocean scenario (not shown), the trends were similar to those of T a , but of weaker magnitude (average of -0.54%/decade). Under this scenario, there is no ice, so the general decrease of the percentage of PAR being transferred by the atmosphere is most likely due to an increase in cloudiness over time. This explanation is consistent with previous studies showing that water vapor content of the atmosphere (Rinke et al., 2019) and cloudiness (Liu et al., 2012) underwent long-term increases related to weaker and less persistent sea ice (Bintanja and Van Der Linden, 2013). For the clear sky scenario (not shown), the trends were of similar magnitude on average than those of the open ocean (average of -0.59%/decade). Under this scenario, there are no clouds, meaning that part of the observed trends in the normalized downwelling irradiance can be ascribed to a reduction in surface albedo. Finally, the magnitude of the trends found for both extreme scenarios was not comparable to the observed trends in T a . This discrepancy can be explained by the contribution of the surface-cloud interactions to the T a trends. We speculate that an indirect consequence of the reduced interactions is causing a change in the spectral quality of the surface incident light over time, as these interactions increase the atmospheric pathlength and thus promotes molecular (Rayleigh) scattering of short wavelengths. This scattering effect could in turn impact the spectral shape of the marine light and therefore have consequences, although modest, on the fraction of PAR absorbed by microalgae.
3.6. Changes in T a compared with possible changes in surface transmittance, T s Areal reduction of sea ice cover and changes in features of the ice pack (thickness, melt ponds, snow cover, ridging, and leads) have been at the forefront of recent Arctic research, notably because they affect the transmittance of the surface, T s , and therefore the amount of light available for microalgae living in sea ice and the water column. That T s is currently increasing in the Arctic is broadly recognized (Nicolaus et al., 2012;Slagstad et al., 2015;Tedesco et al., 2019). These changes in Arctic sea ice are thought to promote more frequent, massive under-ice phytoplankton blooms (Arrigo et al., 2012;Lowry et al., 2014;Assmy et al., 2017) and have been found to trigger earlier and more intense spring phytoplankton blooms (Kahru et al., 2016;Renault et al., 2018). Here, we have assessed that T a decreased over the years for most regions and periods considered in this study. Unfortunately, we could not directly compare the magnitude of this decrease with the increase in T s , because trends in T s at a pan-Arctic scale are not yet available due in part to the difficulty of deriving physical and/or optical properties of snow and sea ice with good accuracy using remote sensing.
To provide a context for the reduction of T a , however, we compiled in situ measurements of T s found in the literature (Figure 8), and the sea surface transmittance for open waters from the model of Gregg and Carder (1990). In summary, for the three main Arctic surface types, namely, snow, bare ice, and melt ponds, T s ranges from 0.2% to 17%, 0.8% to 25%, and 13% to 65%, respectively. Different authors also reported transmittance values for less commonly measured surface types, such as a sea ice ridge, having values ranging between 0.1% and 8.5%, and refrozen leads, having values ranging between 5% (including a snow cover) and 86%. These data are from first and multiyear ice and were collected over regions corresponding to Chukchi and Beaufort Seas, the Canadian Archipelago and Baffin Bay, the Greenland and Norwegian seas, and the Central Arctic Ocean. Additionally, Katlein et al. (2019) published a large collection of spectral values of sea ice transmittance acquired using a remotely operated vehicle, mostly in the Central Arctic Ocean (permanent sea ice). Over our studied period and spectral range of interest, a distribution of their 19,815 T s values is displayed (surface types were not specifically characterized) in Figure 8.
As shown in Figure 8, T s ranges between less than 1% (snow-covered sea ice) and 95% (open waters, leads), with a mean for values from Katlein et al. (2019) below 10%. Across the Arctic, T a has a mean of 60% with extremes of 25% and 75%. Light availability for under-ice phytoplankton is, broadly speaking (referring to the contrasting modes of the data distribution from Katlein et al., 2019, and present study distributions), dominated by the variations in T s . In some places, however, where T s is high, T a could be the primary source of variations of PAR (open waters, refrozen leads, and sea ice with well-developed melt ponds).
During spring, a first drastic increase in T s occurs when the snow starts to melt (the melt onset date generally occurs in early June), and a second drastic increase occurs when the melt ponds spread (a few weeks later in the melt season; Figure 1). Both are key transitions for microalgae (Palmer et al., 2014;Hancke et al., 2018). In the former, T a is generally high (>70%; Figure 6) and likely plays a minor role, but in the latter, T a is lower (45%-60%; Figure 6) and this decrease could partly compensate for the increase in T s .

Conclusion
We gathered a large data set of atmospheric and surface Arctic conditions, which, when applied in radiative transfer simulations, allowed us to assess the variations in the percentage of PAR transferred through the atmosphere (also named normalized downwelling irradiance, T a ) at a pan-Arctic scale between the years 2000 and 2016. Our calculations showed that the maximum PAR available at the surface of the Arctic Ocean generally occurred up to 2 weeks before the summer solstice. In addition, we found that every year, for the permanent ice and transient ice regions, the normalized downwelling irradiance decreased as the melt season progressed, mostly because of the decrease in surface albedo and the ensuing weaker radiative interaction with the overlaying atmosphere. In the case of the permanent open ocean, clouds clearly showed an important impact on the visible light availability throughout the whole solar peak months. In fact, we determined that across the Arctic, if there were no sea ice, E d (PAR, TOA) when propagated through the atmosphere would generally be reduced by 50% due only to the prevalent cloudiness ( Figure 6). This strong impact is dampened by the radiative interaction between the clouds and the sea ice. Finally, we found that on average, the percentage of PAR transferred by the Arctic atmosphere decreased at a statistically significant rate of 2.3%/decade. The spatial domain with the largest negative trends in our normalized downwelling irradiance was the transient ice region. This change over the years probably altered the interaction between the surface and the atmosphere, perhaps modifying the spectral characteristics of marine light.
Most studies of the seasonal evolution of phytoplankton in the Arctic do not identify cloudiness as a determinant factor influencing the annual cycle of light. Although many have emphasized the influence of snow, melt ponds, leads, and ice thickness (Mundy et al., 2005;Arrigo et al., 2012;Assmy et al., 2017;Horvat et al., 2017), much less is known about the variations in the propagation of PAR through the Arctic atmosphere. Here, we found that between 60% and 70% of the top-of-the-atmosphere incident PAR typically reached the Arctic surface, whereas the surface transmittance (T s ) often ranged between 0% and 20%, except for melt ponds which increase T s up to 65%, a value comparable in magnitude to the percentage of PAR being transferred by the atmosphere. Thus, a next step would be to use our framework (Equation 3) to monitor the relation between the atmospheric and surface optically active components (mainly water in diverse states) over a small area at high temporal resolution during the key transition period when E d (PAR, 0 -) increases rapidly due to melt pond formation. Because of the important impact of clouds on light availability, we speculate that on certain occasions, long-lasting clear sky events (which were not discernible with our coarser spatial resolution) could play a role in triggering the phytoplankton spring bloom.
Clouds influence sea ice (Kay et al., 2008;Kay and Gettelman, 2009;Huang et al., 2019;Choi et al., 2020), and sea ice exerts control on clouds (Wang and Key, 2003;Budikova, 2009;Sato et al., 2012;He et al., 2019). This complex feedback interaction is also found when considering their respective influence on the propagation of PAR through the large-scale Arctic atmosphere-ice-ocean system. Further investigation is required to address how, in the warmer Arctic climate, the negative trend in the percentage of PAR transferred through the atmosphere would compensate for the increased surface transmittance. Because both are evolving in opposite directions every season, and also over the years, retrieving T s at the pan-Arctic scale would allow a characterization of the balance for these changes. As a large part of the incident visible light is scattered upward after reaching sea ice, a robust satellite remote sensing method to derive sea ice visible albedo over the Arctic Ocean could help constrain T s .

Data accessibility statement
All data from the model are publicly available at the polar data catalogue (https://www.polardata.ca/pdcsearch/ PDCSearchDOI.jsp?doi_id¼13179).

Supplemental file
The supplemental file for this article can be found as follows: Table S1. Docx