Four Climate Model Intercomparison Project Phase 6 (CMIP6) models were chosen as input to the radiative transfer model (RTM) calculations based on their evaluated ability to reproduce key aspects of high-latitude dynamics and characteristics. The RTM required input from 12 variables (listed below) from each CMIP6 model across two climate scenarios. The selected models include six realizations from the Canadian Earth System Model (CanESM5), four and two realizations from the Max Planck Institute Earth System Model, low and high resolution (MPI-ESM1-LR, MPI-ESM1-HR), respectively, and four realizations from the U.K. Community Earth System Model (UKESM1-0-LL) (Supplementary Table 1). For each combination of model and realization, we extracted and concatenated historical data and future projections for two climate scenarios (SSP2-4.5 and SSP5-8.5) to create continuous time series spanning 1979 to 2100. This process ensured consistency in the data across the entire temporal range, allowing for a comprehensive analysis of climate trends over more than a century (Supplementary Fig. 1). Regarding the choice of ensemble members (Supplementary Table 1), for the CanESM5 model we use both ‘p1’ and ‘p2’, and for UKESM we use ‘f2’ rather than the standard ‘f1’. For all other datasets, we use the standard configuration ‘rXi1p1f1’ where only random factors are varied. For CanESM5, the two physical configurations differ in the remapping method for wind fields (bilinear vs. conservative26), resulting in minor differences, strongest over Antarctica26. Including ensemble members with both ‘p1’ and ‘p2’ physics allow us to explore more of the underlying modelling uncertainty. For the UKESM1-0-LL model, the ‘f2’ configuration is the standard as the ‘f1’ setting is not used66.
To minimize uncertainty from model variability, multiple model responses to anthropogenic forcing were weighted according to model performance and independence. The weights were calculated considering the models’ performance and independence with respect to multiple observational estimates, including climatology, trend, and standard deviation. Performance weights were calculated using the observational records from NOAA Extended Sea Reconstructed SST v5 (ERSSTv567) and the ocean in-situ dataset Coriolis Ocean database for ReAnalysis (CORA5.268). The CORA5.2 dataset includes temperature and salinity in-situ observations both at the surface and at depth while ERSSTv5 is a global monthly analysis of SST data derived from the International Comprehensive Ocean–Atmosphere Dataset. Weights were calculated using the approach described by Lorenz et al. (2018)69 where the model performance is accounted for in the enumerator and in the model independence in the denominator:
$${w}_{i}=\frac{{e}^{\frac{-{D}_{i}^{2}}{{\sigma }_{D}^{2}}}}{1+{\sum }_{j\ne i}^{M}{e}^{\frac{-{S}_{{ij}}^{2}}{{\sigma }_{S}^{2}}}}$$
(1)
In Eq. (1), Di is the distance of model i to observations in space for the three diagnostics (climatological means, trends, standard deviations), \({\sigma }_{D}\) defines how strongly performance is weighted, M the number of models and model realizations, Sij the distance between models i and j, and \({\sigma }_{S}\) defines how strongly model similarity is weighted. The observational ERSSTv5 and CORA5.2 datasets were interpolated to a Cartesian fixed grid of 1×1 longitude-latitude identical to the model grid. Performance metrics were quantified for each diagnostic between modeled and observed sea surface temperature (tos) for the high latitudes (66–80°N) and used to quantify the model weights. Sea surface temperature was chosen as it reflects both direct and indirect metric model skill in seasonally sea ice covered waters. Multiple model realizations of the same model were averaged as part of the weight calculations. Overall, a total of 16 models and realizations were used as forcing of the RTM calculations (Supplementary Table 1). The weighted ensemble averages were calculated for each variable (e.g., PAR, UV, UV-B). The weighted ensemble averages were further weighted for calculations involving area averages. The weighted area average for each latitude-longitude cell was calculated as \(\delta A={R}^{2}\delta \varphi \delta \lambda \cos (\varphi )\) where R is the radius of the Earth, \(\varphi\) is the latitude, and \(\delta \varphi\) and \(\delta \lambda\) are the latitudinal and longitudinal spacing between each cell, respectively.
Calculating spectral irradiance
Irradiance is estimated spectrally using a simple model for the atmospheric radiative transfer of sunlight under clear sky conditions based on specified atmospheric conditions. The spectral band calculations include the wavelengths from 200 to 2700 nm at 10 nm intervals for DNI and DHI, and global horizontal irradiance (GHI) in W m−2 nm−1 falling on a surface horizontal to the surface of the Earth. Calculations follow the ‘simple solar spectral model‘ for direct and diffuse irradiance on horizontal planes at the Earth’s surface implemented as a module in the pvlib70 library. To account for cloudy sky conditions, we apply a ‘cloud opacity factor‘, where we assume that the radiance of a partly cloudy sky can be estimated as a weighted average of the clear sky and overcast sky. The spectral incoming light was initially determined using the Bird Simple Spectral Model (SPECTRL271). Subsequently, we applied monthly average cloud data (representing cloud coverage per grid cell and available as CMIP6 output) to calculate the relative fraction (rho) of diffuse vs direct sunlight using the Campbell-Norman70,72 irradiance equations. This fraction was then applied as described by Ernst et al. (2016)14 to adjust the diffuse and direct spectral light to account for the influence of clouds. The spectral calculations require several input variables. Some were provided directly from CMIP6 model results, such as cloud cover, water vapor content, and ozone thickness. Others, such as albedo, were calculated independently using other CMIP6 outputs. A few variables relied on default parameterizations suggested by the pvlib library70. Below, we describe in detail the individual variables and calculations.
RTM input data
Monthly resolved model outputs from four CMIP6 models and multiple realizations (Supplementary Table 1) were used as input to the light calculations. We calculated the irradiance reaching the ocean surface water (upper 10 cm) under two climate scenarios, SSP2-4.5 and SSP5-8.5, for each combination of CMIP6 model and realization. The RTM, used for the light calculations, required 12 variables extracted from the climate models. The physical and biological model variables included (variable_id: description; Fig. 1, Supplementary Figs. 1–2): tas: air surface temperature, siconc: sea ice concentration; sisnthick: sea ice snow thickness, chl: chlorophyll content, sithick: sea ice thickness, sisnconc: sea ice snow concentration, clt: cloud cover, uas: surface wind velocity east, vas: surface wind velocity north, toz: total ozone column (measurement of total ozone within the atmospheric column), tos: temperature ocean surface, prw: precipitable water (integrated water content of the air). All calculations were done using the CMIP6 datasets made available as Zarr archives by the Pangeo community on Google Cloud73, allowing for calculations without downloading terabytes of data. Prior to calculations, all the required datasets (variables) for each CMIP6 model and ensemble member were interpolated to a Cartesian fixed grid of 1 × 1 longitude-latitude at monthly temporal resolution using a weighted bilinear algorithm74. The RTM was run for each CMIP6 model and realization combination, and scenario independently. The model calculated the light conditions for each grid point (1 × 1 degree) between 60–85°N and 0–360°E for every 4 h (6 timesteps per day) for the 15th of each month between the years 1979–2100.
Ozone data
The ozone fields used in the calculations for the PAR and UV-B were obtained from the CMIP6 input4MIPs data archive75,76. The dataset was the CMIP-recommended ozone forcing for use in CMIP6 climate model simulations that did not represent atmospheric chemistry interactively. The ozone fields consist of a historical simulation (1850–2014) and different shared socio-economic emissions scenarios for the future (2015–2100). For this study, SSP2-4.5 and SSP5-8.5 were used. Both historical and future ozone fields were generated based on simulations of two chemistry-climate models, the US NCAR WACCM-CESM and the Canadian ECCC CMAM models, which participated in the SPARC/IGAC Chemistry-Climate Model Initiative Phase-175,76.
Calculating total column ozone (toc)
Ozone strongly affects the amount of light in the UV spectrum that passes through the atmosphere and enters the ocean. This effect is accounted for by scaling the UV light according to a simple relationship between the thickness of the ozone layer and the amount of UV radiation passing through. Not all climate models include an atmospheric chemistry component and therefore do not provide the total column ozone (toc) as an output variable. However, most models use the same boundary condition and forcing files input required for the CMIP6 various scenario runs (https://esgf-node.llnl.gov/projects/input4mips/), which enables us to use a common ozone dataset across the models in the calculations shown here. The ozone dataset contains ozone volume mixing ratios [mol mol−1] for 1950–2100 for each scenario SSP2-4.5 and SSP5-8.5, which was converted to toz in Dobson Units (DU) at each grid point (x,y,z) and time (t) for use in spectral irradiance calculations. The conversion was done using the equation:
$${toz}=10 \cdot \frac{\left(R{T}_{0}\right)}{\left({g}_{0}{P}_{0}\right)} \cdot \left({\sum }_{i=1}^{N-1}0.5({VMR}\left(i\right)+{VMR}(i+1))\left(\right.p\left(i\right)-p(i+1)\right)$$
(2)
where VMR is the mixing ratio (ppm), N is the number of vertical levels of the air column, R = 287.3 is the specific gas constant for air (J kg−1 K−1), T0 = 273.15 is temperature (K), P0 = 1.01325e5 is the standard pressure at surface (Pa), g0 = 9.80665 is the global average gravity at the surface (m s−2), Na = 6.0220e23 is Avogadro´s number, p is pressure in hPA. Calculated values for the period 1/1/1950 to 12/31/2099 ranged over 218.9–614.7 DU with a mean of 332.1 DU for SSP5-8.5. For scenario SSP2-4.5, values ranged over 218.9–558.8 DU with a mean of 326.6 DU.
Calculating ocean surface, sea ice, and snow albedosOcean surface albedo
The albedo over the open ocean is usually a constant equal to 0.06 in climate models and rarely are more complex approaches used. A recent publication13 proposes a new approach for next-generation climate models that quantifies the effects of the solar zenith angle, ocean waves, and chlorophyll content on ocean surface albedo (OSA), which can help reduce the uncertainty in climate sensitivity to the flow of radiative energy. The proposed OSA scheme, which was implemented in the RTM presented here, calculates the various contributions spectrally from the ocean surface on both direct and diffuse shortwave radiation, providing a more realistic approximation of the reflected shortwave radiation. The complex implementation is described in Séférian et al. (2018)13 and not repeated here.
Sea ice and snow albedo
Sea ice and snow are very efficient mediums for reflecting shortwave radiation and one of the most important factors influencing the Arctic energy budget. However, the efficiency can vary with sea ice thickness, snow crystal structure, and melting ponds, or the purity of ice such as soot particles originating from anthropogenic dust, black carbon from coal combustion, or volcanoes. Here we consider sea ice and snow as pure as we do not have information on the geospatial distribution of any properties that could decrease albedo in the future except purely physical changes. The parameterization of the annual average albedo of thick (hice > 0.5 m) ice was 0.52 and albedo of snow-covered ice was 0.65. These values are based on the Community Climate Model System Version 3 (CCSM377,78) and the Community Ice Code (CICE12) component of the Community Earth System Model and modified to reflect the seasonal changes in sea and snow conditions. The values used for albedo of snow and ice reflect average annual conditions that include seasonal deterioration of sea ice and snow surface conditions away from pure conditions (ice = 0.73 and snow = 0.96). These values are comparable with recent observed values where average albedo of Arctic sea ice was found to be 0.8 in April–May and decreased to 0.4 between June and August41,79. When ice thickness is less than 0.5 m and the ice is not melting, the dry albedo can be defined as \(a({dry})={a}_{o}\left(1-{f}_{h}\right)+0.52{f}_{h}\) where \({a}_{o}\) is the open ocean albedo and \({f}_{h}\) is an asymptotic function defined by \({f}_{h}=\min \left(\frac{{tan }^{-1}\left({c}_{{fh}}h\right)}{{tan }^{-1}\left({c}_{{fh}}0.5\right)},1.0\right)\) where \({c}_{{fh}}=5.0\) and h is the ice thickness80. To represent impact from melt ponds on the albedo, we used the near-surface air temperature (tas) and calculated the wet albedo as \(a\left({wet}\right)=a\left({dry}\right)-{0.075f}_{T}\) where \({f}_{T}=\min \left({T}_{{air}}-1.0,\,0.0\right)\;{{\rm{and}}}\,{T}_{{air}} > -{1.0\,}^{\circ }{{\rm{C}}}\). The algorithms for sea ice and snow can also be found as part of the UCAR Community Earth System Model equivalent to CESM-CICE 5.0 (https://github.com/CICE-Consortium).
Attenuation from snow, ice, and chlorophyll-a
Earlier observations have shown that the attenuation coefficient through snow can vary considerably, ranging from 4.3 to 40 m−1 18, and ice modeling has reflected this. For example, the CCSM3 model assumed that no shortwave radiation penetrates the snow80, while the ROMS sea ice module used 20 m−1 by default81. The latter value was used in this study as it represented a value in the middle of the range of observed values18, and was also used in Castellani et al. (2022)4. Attenuation of snow on top of sea ice was estimated as an exponential decaying function82 of snow thickness (snthick, \({h}_{{snow}}\)) using a fixed coefficient of \({k}_{{snow}}=20\) [m−1]: \({I}_{{ice}}={I}_{{sfc}}{e}^{-{k}_{{snow}}{h}_{{snow}}}\) for all wavelengths.The spectral absorption coefficients of sea ice describing how light is attenuated as it propagates through the ice was estimated for the wavelength range 200–1000 nm by combining previously published observations82 that were interpolated to a fixed wavelength of 10 nm. Total attenuation of light from the surface (top of ice, but underneath snow if present) to the underside of the ice was calculated for each wavelength (\(k(\lambda )\)) as an exponential function82 of the thickness of sea ice \({I}_{{sfc}}\left(\lambda \right)={I}_{{ocn}}\left(\lambda \right){e}^{-k\left(\lambda \right){h}_{{ice}}}\). For these calculations, we assumed pure ice, and we did not account for any black carbon (e.g., anthropogenic sources of soot from coal combustion) that would affect the ice optical properties of absorption and reflection. We also do not account for attenuation caused by sea ice algae83, as the CMIP6 models do not provide that information. We extracted spectral absorption coefficients of phytoplankton (chlorophyll-a) from Table 3 of Matsuoka et al. (2007)19 and used them to quantify the exponential decay of light with depth due to chlorophyll-a. The chlorophyll in the surface waters results in a further attenuation factor that is independent of the attenuation caused by other seawater components (the water itself, cDOM, other particles). We have assumed that the attenuation from chlorophyll, over a layer (0–5 m), should be independent of the depth distribution within that layer84.
Comparison with ERA5 and CMIP6 shortwave radiation
When comparing the monthly averaged global horizontal incoming (GHI) shortwave radiation in our model with ERA585 for the period 1979–2020, we found a small, consistent, spatially homogenous bias. This bias can be explained in part by the coarse resolution in our model inputs (1 × 1 degree longitude-latitude) compared to the high resolution of ERA5, and that our RTM model is simpler in capturing the most essential elements, or processes, required to quantify shortwave radiation. Our RTM also lacks several important atmospheric components, such as detailed cloud layering, cloud thickness, and reflection from land, among others. There is also a difference in the land/ocean mask which makes a difference in albedo between the ocean and land between the RTM and ERA5 (not shown). We adjusted for these systematic errors by multiplying GHI by a factor of 1.17. After bias-correcting, the modeled surface incoming shortwave radiation (GHI, 200–2700 nm) had a correlation of r = 0.998 (p < 0.05, n = 480) and RMSE = 2.8 W m−2 with the monthly mean surface shortwave radiation from ERA5 global reanalysis for the period 1979–2020 (Supplementary Fig. 10).
Next, we compared RTM outputs of absorbed shortwave radiation (incoming shortwave minus effect of albedo) at the surface with those estimated by the individual CMIP6 models (CMIP6 table_ids: rsds, rsus) and found that the RTM skillfully replicates the values, ranges, and the dynamical variability between 1979 and 2100. For the Barents Sea LME, the average temporal correlation was r = 0.95 (p < 0.05, n = 60) and RMSE 4.6 W m−2 across 16 CMIP6 model datasets between 1979 and 2100. For the Northern Bering and Chukchi Seas, the correlation was r = 0.91 (p < 0.05, n = 60) and RMSE = 6.5 W m−2 (Supplementary Fig. 12).
The RTM closely aligns with individual CMIP6 models in calculating incoming GHI light and effectively captures the trend and dynamic variability from 1979 to 2100, demonstrating a high correlation (Supplementary Table 2). Particularly, the RTM replicated the estimated light of the CanESM5 and MPI-ESM2-HR models closely, while the UKESM1-0-LL had a higher bias (Supplementary Fig. 13). The RTM effectively replicated trends and relative values across two marine ecosystems with significant seasonal fluctuations in sea ice extent and light conditions, thereby confirming the reliability of our projections of relative changes in light levels.
Model sensitivity
Some of the RTM model parameters vary considerably when considering the literature, such as the choice of snow attenuation coefficient, which is reported between 4.3 and 40 m−1 18. To understand the importance of some key choices we made for the RTM, and how these affected the final outputs, we performed a sensitivity test. We ran the model for 10 years (1/1/1979 to 1/1/1989) with either a feature turned on or off or a change in parameter value. Features that were turned on or off included the effect of melt ponds (impacts albedo), wind (impacts surface roughness and albedo), chlorophyll (impacts albedo and attenuation), sea ice concentration (impacts attenuation and albedo (not shown in Supplementary Fig. 11)), and the effect of using a snow attenuation coefficient of 5.9 m−1 versus the default 20 m−1. We also estimated the effect of using the more complex ocean surface albedo (OSA) scheme versus a default value of 0.06 for ocean albedo. We analyzed the average percent difference between the features turned on or off for each individual model and model member in outputs from the RTM over a 10-year simulation (Supplementary Fig. 11). Overall, the mean differences with or without a model feature changed the amount of PAR reaching the water column by <14% for both the Barents and Northern Bering and Chukchi Seas. Change in snow attenuation parameter increased PAR entering the surface waters by an average of 8% (Supplementary Fig. 11) for both LMEs. Removing the effect of chlorophyll (0.05%), and surface wind (1.0–1.16%) also increased the amount of light reaching the surface of the water column by reducing albedo and attenuation. Removing melt ponds reduced the amount of PAR entering the ocean by 0.05% for the Barents Sea and 0.13% for the Northern Bering and Chukchi Sea. The effects from individual features on UV-B were comparable to the effects on PAR, although the relative effect was smaller with changes <7.5%. The effect of implementing a spectral method to calculate ocean surface albedo, OSA, increased the albedo relative to a fixed ocean value of 0.06 with an average impact of reducing PAR by 1.3–1.6% (Supplementary Fig. 11). Generally, the sensitivity tests suggested that changes in our approach or choice of parameter can contribute to changes of up to 15% impact on the annual average amount of PAR that reaches the water column.
Model limitations
Our analysis only considers the light levels just below the sea surface and does not account for changes in attenuation within the water column, e.g., due to changing concentrations of pigments, colored dissolved organic matter (cDOM), and particulates in part driven by changing terrestrial input. Future work may propagate irradiance down the water column using CMIP model output, which partially covers the previously mentioned variables. However, noting that much of the changes in pigment concentration may occur within subsurface maxima and assuming that most of the newly exposed sea area is not strongly impacted by coastal inputs, it seems likely that the first-order impacts on surface-layer irradiance are already captured by changes in atmospheric attenuation and ice cover/sea-surface albedo as calculated herein. Still, under-ice blooms during spring could impact the attenuation of light reaching the water column, which the model does currently not account for. The effect of clouds, snow, and ice dynamics are all calculated as one dimensional, which reduces the complexity of the RTM, but also means that we are simplifying multi-layered dynamics. A future version of the RTM may include more realistic layering to better resolve such features.
Biological model
To model the survival of eggs and growth of juveniles of different species, we used functional relationships with temperature as observed in laboratory experiments34,43,86. Temperature for the winter months (February, March, April) were used as input to the functional relationship between egg survival and temperature for Atlantic cod, polar cod, and walleye pollock. Summer (May, June, July) temperatures were used to derive the juvenile growth rates for each species. Calculations of either egg survival or juvenile growth were performed for each grid point within the LMEs and then averaged to provide a timeseries with regional uncertainty. The temperatures and values that maximize egg survival and juvenile growth for each species can be seen in Table 1 (optimal temperature). This simple biological model did not consider direct impacts of light on feeding, only the indirect effects of changing light conditions on temperature in the water column.
Ethics and Inclusion
Our authorship team included one researcher who was based in the Barents Sea region throughout the study design, implementation, data analysis, and partly during the authorship stages. Our results are designed to inform regional rather than local-scale activities, Regional research relevant to this study has been appropriately cited.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.