Satellite-based evidence of regional and seasonal Arctic cooling by brighter and wetter clouds

. Using two decades of satellite-based measurements of reﬂectance of solar radiation at the top-of-atmosphere and a complementary record of cloud properties, it is concluded that the loss of Arctic brightness due to sea ice retreat is not compensated by a pan-Arctic increase in cloudiness, but rather by a systematic change in the thermodynamic phase of cloud and a resultant effect on cloud reﬂectance. Liquid water content of the clouds has increased resulting in positive trends in susceptible cloud properties. Consequently, a cooling trend by clouds is superimposed on top of the pan-Arctic ampliﬁed 5 warming, induced by the anthropogenic release of greenhouse gases, the ice albedo feedback and related effects. Except above the permanent and marginal sea ice zone around the Arctic circle, the rate of surface cooling by clouds has increased, both in spring ( − 32% in total radiative forcing for the whole Arctic) and in summer ( − 14%). The magnitude of this effect depends on both the underlying surface type and changes in the regional Arctic climate. the other Arctic show a negative total CRF, − 51.1 W m − 2 over Greenland Sea and − 42.2 W m − 2 Barents in AMJ to the − 39.2 W m − over those regions inﬂuenced by the climate of low latitudes (Bafﬁn Bay, Greenland and Barents and − 30.7 W m − 2 and − 22.4 W m − over the Hudson Bay and Kara Sea in JAS, respectively. that implies an increasing amount of supercooled liquid cloud droplets. The higher reﬂectance of clouds results in a more negative radiative forcing at the surface, thereby locally dampening Arctic Ampliﬁcation. However, cooling by clouds implies the strengthening of the meridional temperature gradient. We expect this will lead to increase the inﬂow of warmer and moister air masses from the lower latitudes into the Arctic climate. This may then either further decrease Arctic Ampliﬁcation by generating more supercooled liquid water cloud, or possibly enhance Arctic Ampliﬁcation by the increase input of warmer air, or some combination of the two. Future model projections of the Arctic climate must take into account these effects to accurately predict the impact of anthropogenic emissions of greenhouse gases and short-lived climate pollutants.

the the perennial sea ice edge shelves to CFC and downwelling LW during spring months.
In addition to the warming from the release of greenhouse gases, the dynamics of air masses and physical properties of clouds may contribute to the tropospheric thermal emission. CFC aside, the optical properties of clouds, such as the optical thickness (COT, τ ) and effective radius r eff of droplets/crystals, and liquid/ice water content (LWP/IWP), regulate downwelling LW. Model projections of the century ahead show that Arctic clouds during summer are rather unaffected by sea ice variability 80 but their response to sea ice loss is to become optically thicker, have higher LWP, and more frequent low-level liquid phase within the boundary layer (Morrison et al., 2019). Hence, the changes in τ and thermodynamic phase of clouds enhance or suppress cloud radiative forcing (CRF) at the surface. This behavior has been identified in the continuous surface measurements above the Beaufort and Chuckchi Seas (Shupe and Intrieri, 2004), at Ny-Ålesund, Svalbard (Ebell et al., 2019) and in the data products retrieved from AVHRR (Francis and Hunter, 2006). 85 From the above review of current knowledge about the changing conditions in the Arctic, we consider that investigations of the R TOA λ and the cloud properties over the past two decades will provide insight into the evolution of the Arctic Amplification.We have prepared a consolidated R TOA λ data set from 1995 to 2018. A key set of satellite sensors record backscattered radiation in the solar portion of the spectrum. Consequently, this study focuses on the months between April and September. The Arctic seasons considered are spring, defined for our purposes as April May June (AMJ) and summer, July August 90 September (JAS) (see App. 2.1). The investigation of R TOA λ trends involved analysis of the last twenty years of cloud properties from the observations of AVHRR, retrieved with the latest algorithms ) (see App. 2.2), which supersede older popular data sets, for which specific errors have been found (Zygmuntowska et al., 2012). We build on the heritage of the earlier studies describing the Arctic state and extend the trend analyses limited previously to 1982-1999(Wang and Key, 2005b. 95 The objectives of this paper are threefold. Firstly, we provide evidence that space-borne measured, spectral R TOA λ is a valuable indicator of the changing atmospheric composition and surface properties of the Arctic. Secondly, we determine R TOA λ trends at regional and seasonal scales and reveal unexpected patterns of behavior. Thirdly, we attribute the trends in R TOA λ to changes in the thermodynamic phase of clouds. Lastly, we quantify the average radiative forcing by clouds and its changes. We relate the latter to the changes in the physical properties of clouds in response to climate change.

2 Data and Methods
The study of the Arctic geophysical domains by remote sensing requires sensors with broad spectral coverage at an adequate spectral resolution to be able to separate different spectral features of gases, surfaces, liquid water and ice/snow. We define the spectral reflectance measured at TOA -R TOA λ -as 105 where I λ is the Earthshine, i.e. the upwelling scalar radiance measured at TOA (units of photons × s −1 cm −2 nm −1 sr −1 ), E 0 λ the unpolarized downwelling solar irradiance (photons × s −1 cm −2 nm −1 ) and θ 0 the solar zenith angle in degrees.
Relevant parameters for the R TOA λ analysis are shown in Fig. 1. The left y-axis of Fig. 1-a shows I λ , E 0 λ for a GOME measurement above the Kara Sea, while the right y-axis shows modeled R TOA λ , in satellite perspective, representing the TOA signal for typical Arctic geophysical conditions. Fig. 1-b shows the wavelength dependence, at GOME spectral resolution, of the spectral 110 reflectance for different surface types. The almost flat Earthshine between 450-800 nm reveals the presence of a cloud deck or snow surface in the satellite field of view. Ten wavelength bands of spectral width 5-10 nm have been selected satisfying the following requirements: they are chosen to be similar to those of sensors' channels used in the literature for comparative purposes. Their coverage from the UV to the NIR shall provide differential sensitivity to individual components of the Arctic atmosphere-surface. Additionally, they should exclude strong absorption by atmospheric trace gases to avoid misinterpretation 115 of the observed behavior. Two exceptions are the broadband O 3 Chappuis band (525-675 nm) and the narrow O 2 A-band (centred at 760 nm). The former, even if smoothed at 5-10 nm resolution, has information on the total column of ozone and on the structure of the upper troposphere and lower stratosphere. Well-mixed gas, such as oxygen, provide valuable diagnostics about the depth of the atmospheric column, as seen from space. The A-band is used to assess the surface topography in a cloudfree atmosphere (van Diedenhoven et al., 2005) and altitude, geometrical and optical depth of clouds over dark (Rozanov and 120 Kokhanovsky, 2004;Lelli et al., 2012Lelli et al., , 2014 and bright (Schlundt et al., 2013) surfaces.

Reflectance data
To detect changes on daily, monthly, seasonal and decadal scales several measurements per day of adequate spatial resolution must be made over several decades. The polar-orbiting spectrometer suite comprising GOME, SCIAMACHY and GOME-2 Figure 1. Plots of the solar irradiance, the radiance of a cloud (Earthshine) and reflectances at top (TOA) and bottom (BOA) of the atmosphere as a function of wavelength from 280 nm to 800 nm. The cloud radiance was observed by GOME on May 15, 2001 over Kara Sea (80.53 • N, 75.99 • E). Modelled R TOA λ (nadir, solar zenith 40 • ) display a water cloud, placed at 3 km and optically dense 30, above sea water and snow, with a cloud-free sea ice, snow and melt pond spectrum. The lower panel shows the black sky hemispherical reflectance at the ground of relevant Arctic surface components. Chlorophyll absorption is taken from Clementson and Wojtasiewicz (2019) and plotted for a May 2016 concentration of 12 mg m −3 observed in the Bering Sea (Frey et al., 2018). Arctic shrub and coarse snow data are taken from the ECOSTRESS and ASTER spectral libraries (Meerdink et al., 2019;Baldridge et al., 2009). Melt pond and sea ice albedos are from Istomina et al. (2013). (Tab. A2 for their specifications) is an optimal choice given the individual length of the time series and their high spectral 125 resolution. Description of GOME can be found in Burrows et al. (1999), while SCIAMACHY and GOME-2 are respectively provided in Burrows et al. (1995) and Munro et al. (2016). The detailed steps to harmonize R TOA λ measured by sensors of different technical specifications are given in the App. A.
Sensors measuring solar radiation scattered back to the top of the atmosphere do not measure at night. Measurements in thermal infrared (λ 4 µm) are required to record the thermal emission from the surface and the atmosphere. In practice, 130 this situation, coupled with different sensor swath widths, limits the use of R TOA λ measurements up to a common north parallel for those months that guarantee adequate sampling. This effect is illustrated plotting the pan-Arctic annual cycle of R TOA λ in Fig. 2. At the three wavelengths 510, 560, and 760 nm, the seasonality shows that summer months have lower R TOA λ and higher otherwise. The darkening of the Arctic can also be seen by comparing the years at the beginning of the record, 1996, Figure 2. Annual cycle of spectral R TOA at three wavelengths (λ = 510, 560, 760 nm) for the full record from 1996 to 2018. All sets exhibit the demarcation between months of steep (Apr-May-Jun) and flat gradient of R TOA (Jul-Aug-Sep). This shift leads by one month melt onset (6 Jun), followed by sea ice opening, breakup, minimum (16 Jul -Sep inclusive), and freeze onset (4 Oct) as observed with satellite brightness temperatures (Smith et al., 2020). On the rightmost panel the terminator location of the three sensors with the 85 • N (grey line) common threshold used for monthly R TOA aggregation.
with the most recent ones. However, this behavior occurs between April and September. These are the months when the the R TOA λ curves show that recent years are brighter than those at the beginning of the time series.
From Fig. 2 we identify two distinct behaviors of R TOA . One of steepest decrease, from April to June, and one relatively flat, between July and September. In contrast to the common definition of the climatological seasons, we group April, May, 140 June (AMJ) as Arctic spring and July, August, September (JAS) as Arctic summer. This distinction is explained by the sensors' measurement strategy and by the time-dependent physical processes leading to the transition between high-to-low Arctic reflectance in June to the minimum sea ice extent in September. The changes in surface reflectance between April and May are attributed to snow cover changes and those in June to sea ice changes (Smith et al., 2020). Over water, the timing of such transitions increasingly approaches the summer solstice, which is the day of strongest solar insolation, while it moves further 145 away from it over land (Letterly et al., 2018). It is reasonable to regard this day as a more natural demarcation point between Arctic spring and summer.

Cloud and flux data
In our study, the R TOA λ data is complemented by a record of cloud properties and fluxes at TOA and BOA. These are inferred from the afternoon orbit of AVHRR sensors onboard the POES missions. The primary reason for choosing these records is 150 the abundance of studies using these data in the Arctic. This has the required coherent radiometric calibration before the implementation of individual downstream methodology to assess changes across the Arctic. The cloud and flux records, version 3, are presented by Stengel et al. (2020). The application of the cloud algorithm to MODIS measurements, mimicking AVHRR channels, has shown that the retrieval scheme is well aligned with the reference standards of CloudSat and CALIPSO data for CFC, CTH, COT and liquid thermodynamic phase. While agreeing on the sorting of cloud tops between water and ice phases, 155 higher IWP variability is found as compared to that in the reference DARDAR cloud data products (Delanoë and Hogan, 2010). Aggregated IWP histograms do not differ (Stengel et al., 2015). Version 3 has improved version 2 in terms of precision, accuracy and stability (Stengel et al., 2017). Even more relevant to our purpose is the scheme adopted to calculate cloud properties and broadband fluxes. CFC and optical property retrieval uses a neural network trained on CALIOP data to account for the extent of the underlying bright Arctic surface. CTH has been corrected with CALIOP profiles to account for the depth of 160 light penetration inside a cloud. The retrievals of CTH from all infrared thermal channels are influenced by this effect and yield a radiative cloud top height, lower than the physical cloud top. Broadband fluxes are not derived incorporating reanalysis data but the retrieved cloud properties instead. This is important, because the changes in fluxes are coherently related to changes in the cloud properties. Given the standard notation (all = all-sky, clr = clear-sky, + = upwelling and − = downwelling fluxes F ), the average long-term relative bias of AVHRR-derived fluxes against CERES ranges from +2.  . Misclassified cloudy scenes especially over dynamically bright surfaces (i.e. marginal and fractional sea ice zones) impact the calculation of broadband fluxes. This has been already noted in first studies comparing ERBE and AVHRR cloud radiative forcing derived with different scene classification schemes (Li and Leighton, 1991). The conversion of directional radiance, 170 measured at TOA, to irradiance requires the knowledge of the angular light redistribution function of the surface and atmospheric components. If this is not accurately assumed, the irradiance and F +/−,clr SW above reflecting surfaces cannot optimally be calculated. Using the same data of our study, it has been found a low sensitivity of trends in cloud radiative forcing to the biases in cloud properties over surfaces of changing brightness (App. d in Philipp et al., 2020, p. 7499). This confirms the suitability of cloud properties retrieved from AVHRR measurements, classified with active sensor data, for Arctic trend studies.

Results
The R TOA λ time series, measured by GOME, SCIAMACHY and GOME-2A over the Arctic region (60-85 • N), anomalies, trends and significance were harmonized (for more details see App. A and App. B). They are shown for wavelengths 510, 560 and 620 nm in Figs. 3 and 4. The R TOA λ retrieved from the sensors MERIS (on Envisat) and GOME-2B (on MetOp-B) confirm that the correction scheme is successful. A consistent and consolidated data set results from the measurements of the three 180 instruments. Any errors are minimized, when sunlight availability across the Arctic provides full coverage for the sensors' swath at highest latitudes. Seasonality is the dominant feature of Fig. 3. Maximum R TOA λ occurs in early AMJ when the Polar day results in the Arctic being fully illuminated and the ice extent is close to its maximum. Analogously, minimum R TOA λ occurs from August to September when the days are shortening and sea ice coverage is at its minimum. The observed seasonal cycle of R TOA λ agrees with that observed by models as do the observations of sea ice extent over the Arctic (Holland et al., 2008).  Morrison et al. (2019) state that no significant relationship between CFC patterns and sea-ice loss is observed during summer but some is identified in autumn months. Such changes are not observable in the pan-Arctic R TOA anomalies. Rather, the reduction in reflectance is small and not attributable to a specific season. As a consequence, we need to ask whether the loss of reflectance associated with sea ice reduction is compensated by increasing CFC or brighter clouds and which processes lead to the small pan-Arctic R TOA trends.

195
To answer these questions in the following, we map R TOA λ in the Arctic, gridded at 1 • × 1.5 • latitude and longitude. Fig. 5 shows the spatially resolved R TOA λ trends for λ = 510, 560, 620 nm over the Arctic region for AMJ and JAS. Similarly, Fig. 6 shows trends for the analyzed wavelengths for the 12 Arctic regions, that are defined using the geographical subdivision proposed by Serreze and Barry (2014) and Wang and Key (2005a)   There are marked regional differences. Those that are statistically significant (at 95% confidence level) are shown with red crosses. For AMJ a significant negative trend over the Barents Sea is compensated at all four wavelength bands, by a positive Although not of the same magnitude, almost all regions show a reflectance change at λ 760 nm. Oxygen absorption is modulated primarily by CTH and, to a lesser extent, by CFC and optical properties such as CA and COT. In this context, which indicate a clear change in occurrence of clouds or one of their scattering properties. Knowing that R TOA is influenced by 220 scattering and absorption in the atmosphere (Sledd and L'Ecuyer, 2019;Donohoe and Battisti, 2011) and that the atmospheric R TOA can be additionally partitioned into cloud, aerosol and gas contributions, this prompted us to examine changes in those cloud properties which directly influence the spectral R TOA trends.

Cloud properties
The globally-validated and consolidated cloud record  has first been analyzed across the Arctic (60 • -225 85 • N). The top plot of Fig. 7 shows time series of CFC and CTH. Both parameters show small, statistically insignificant, trends over the last 20 years. CFC has slightly increased by about 0.001 (+0.14%) decade −1 while cloud tops are lower by ≈ 6 m (−0.14%) decade −1 . This finding obviously excludes an explanation being that reflectance loss at visible wavelengths, due to shrinking sea ice extent, is offset by more CFC or that the loss of CFC reveals more bright underlying surfaces. However, the middle plot of Fig. 7 shows that over two decades the COT temporal trend of liquid clouds has opposite sign of that of cold Similar to the approach we used for the R TOA λ trends regionally and qualitatively, we map cloud parameters in the bottom panel of Fig. 7, adding also the albedo of clouds at λ = 600 nm. CFC trends are regionally partitioned and are seen to increase in the range 5-20% where the greatest sea ice losses are observed. This occurs during AMJ and less extensively in JAS.

235
Examples of this behavior are found in the Barents, the Kara and Laptev Seas. On the contrary, large areas of statistically significant decrease in the range 2.5-10% are homogeneously observed across land masses circling the inner polar belt. This includes Greenland and the Atlantic corridor, confirming past results (Hofer et al., 2017). More pronounced trends of the 11 https://doi.org/10.5194/acp-2022-28 Preprint. Discussion started: 24 January 2022 c Author(s) 2022. CC BY 4.0 License. different cloud parameters, irrespective of their sign, occur in AMJ rather than in JAS. The Hudson Bay is one of the few regions experiencing a seasonal trend reversal. The AMJ period is characterized by less cloudiness (-5%), whereas the JAS 240 period exhibits an increase of the order of almost 10% over the last two decades. The resemblance to the trend reversal of all R TOA channels (Fig. 6) indicates that CFC changes primarily modulate R TOA λ over the Hudson Bay. This is inferred from the absence of change of trend sign of those cloud parameters that influence the reflectance in the solar spectrum, such as COT of water and ice clouds (third and fourth column in Fig. 7). decades CTH in these regions has decreased by 10% on average. In JAS, however, CTH increases significantly from the Fram strait, throughout the Barents and Laptev Seas, closer to the Pole, and western Siberia, with a less pronounced positive trend for Greenland (+4%) and a more pronounced one for the Hudson Bay (+8%), albeit non significant. Total COT is split into liquid and solid cloud phases. The geographic distribution of the trends in Fig. 7 provides insight into which areas are responsible for 250 the positive pan-Arctic trend in τ of liquid clouds (τ -water) and for the negative trend for ice clouds (τ -ice). τ -water increases across the whole Arctic in AMJ except over the Atlantic sector and the southern part of the Hudson Bay. The positive trend is maintained over north Greenland, Canadian Archipelago, North Pole and part of the Eurasian continent also during JAS.
A positive trend of τ -water is correlated with a trend of opposite sign for τ -ice: this holds for all regions of permanent and sector show a different behavior: there is a 34% increase in τ -water during AMJ and 22% increase in JAS. Notwithstanding the increase over localized areas (e.g. north Greenland), mean τ -ice over the Arctic regions remains nearly unchanged in both seasons. The liquid phase of clouds does not increase across the Fram Strait, whereas the ice phase decreases by roughly 20% in both AMJ and JAS periods. Finally, the Atlantic sector (the Greenland and the Norwegian Seas) show decreases in the COT for both the liquid and solid cloud phase during AMJ and JAS.

260
The rightmost polar plots of Fig. 7 show seasonal trends in cloud albedo (CA), for which a marked change of the spatial rather than temporal scale is observed. To a certain extent, the CA trends are geographically correlated with those of CFC and τ -water.
Individual regions are grouped in a similar manner to the R TOA polar plots: comparable distribution of CA are found over the most eastern and most western Arctic Seas (Beaufort and Chuckchi, East Siberian, Laptev, and Kara Seas). Positive trends are almost invariably distributed over water masses, the Canadian Archipelago and the northern part of Greenland, irrespective of 265 the season. In contrast, clouds become darker at lower latitudes, southern Greenland and the Atlantic sector. Over the Siberian land masses this is not observed, and CA changes in the region, besides being small, are attributed to a competition between changes in CFC and τ -water. The loss of albedo due to cloud dissipation is compensated by the increment in albedo through increased τ -water.
To facilitate a quantitative seasonal comparison between the Arctic sectors, Fig. 8 shows the trends and the standard error 270 (i.e. 2-σ standard deviation, see App. 2.2) of five cloud properties (CFC, CTH, τ of water and ice phase, CA) together with trend of liquid (LWP) and ice water path (IWP), from the same cloud record . Changes in R TOA λ depend in the first place on changes in cloudiness and τ (irrespective of the phase), which in turn is a function of LWP, droplet/crystal effective radius (r eff ) and air density ρ (i.e. τ = 3/2 × LWP/ρ r eff ). The sign of LWP and IWP trends confirm the τ trends. We infer that τ -water has increased as a result of the positive change of LWP and/or a concurrent systematic pan-Arctic decreasing 275 trend of r eff (see Fig. C1).

Cloud radiative forcing
We compute the net radiative forcing,F cld or CRF, due only to clouds. The flux data are taken from the cloud record ) (see App. 2.2). The difference between F − and F + is the net radiationF and CRF =F all −F clr . The multi-year mean and trends of SW, LW and total CRF at the surface are plotted in Fig. 9.

280
At pan-Arctic scale clouds exert a negative SW radiative forcing of −58.7 and −63.8 W m −2 in AMJ and JAS, respectively.
In the same seasons, the LW component amounts to +46.9 and +46.1 W m −2 , and the multi-year mean of total CRF is −11.8 and −17.6 W m −2 . However, CRF is seasonally and regionally partitioned: clouds' total radiative forcing at the surface is positive over bright areas as a result of LW offsetting SW. For instance, total CRF over Greenland is +14.9 and +23.5 W m −2 , which corresponds to the Arctic sectors over which clouds reflect the least in absolute SW (−19.8   From the CRF trends of the last two decades (Fig. 9), clouds over the perennial sea ice zone are increasingly cooling TOA (see Fig. D1) and BOA alike, while being neutral to positive over the Atlantic corridor and land masses at low latitudes. In AMJ

Discussion and conclusions
In the last two decades, the set of analyzed parameters provides a coherent geophysical picture: the Arctic R TOA λ has declined.
However, this decline is less than that expected as a result of the loss of sea ice. We attribute the reason for this decreasing  To some extent, Wang and Key (2005b) anticipate the results of our work. The downward trend in broadband albedo of −1.40% decade −1 between 1985 -1999 is confirmed by our negative all-sky R TOA λ trends, implying a sustained sea ice loss after 2000 and general darkening of the Arctic surface. However, the regional patterns match neither our results nor most recent knowledge (Hofer et al., 2017). The annual increase of 0.6% in CFC over Canadian Archipelago, Chuckchi Sea and Siberia and, in JAS, over Greenland reported in Wang and Key (2005b) is probably explained by the limited length of the analysed 325 record. Trends in CFC over Greenland, for instance, level out before 1995 but turn strongly negative afterward, contributing to a significant loss of the ice shield mass (Hofer et al., 2017). This explains the nonexistent clouds' τ trend in Wang and Key (2005b), which is in contrast to the significant moistening across most of the Arctic of Figs. 6, 7 and 8.
Greenland has a unique behavior: R TOA λ trends at all wavelengths are positive, irrespective of the season (Fig. 6). The AMJ R TOA λ trends, up to 5%, are even larger than those for JAS. This result is particularly surprising, given the insignificant CFC 330 trend (Fig. 7,8), thus not contributing to an increase of the overall reflectance. Similar behavior is found in the Hudson Bay and Canadian Archipelago, which show an increase in reflectance, in contrast to a general darkening of the Arctic. The mechanism by which these regions increase R TOA λ lies in the link between LWP and CA, through τ -water. In fact, τ -water changes sustain the correlated R TOA λ changes because of the non-linear relationship of CA to τ -water via LWP. It follows that a R TOA λ loss is overcompensated by wetter clouds in the northern sector and by increased snowfall in the southern part of the Greenland 335 continent. Cloud LWP has increased by 28-30% over Greenland and by 14-16% over the Hudson Bay and the Canadian Archipelago, having positive τ -water trends of 30%, 14% and 22%, respectively. Notably, the seasonal behavior of τ -water, increasing over Greenland, is not associated with CFC loss and a positive CRF change in the last 20 years. In contrast, cloud dissipation, increased by anticyclonic activity and concurrent temperature inversion strengths, is responsible for enhanced insolation at the ground and melting (Hofer et al., 2017). JAS R TOA 560 changes over the Hudson Bay are exceptional. They are 340 correlated with a 9% increase in τ -water and minimal CRF changes. This area shows one of the largest CFC increases during JAS months (Fig. 8), also corroborated by similar significant changes in AMJ and JAS observed in the reanalysis data (Fazel-Rastgar, 2020). The total CRF is −30.7 W m −2 while +2.7 W m −2 during AMJ. CRF trends point to a cloud cooling of the Hudson Bay at a rate of −2.9 (AMJ) and −1.3 (JAS) W m −2 over the last two decades.
Cloud forcing at the surface depends on cloud property changes. The behavior is summarized in the seasonal and regional 345 charts of Fig. 10, in which mean value and trend of SW, LW and total CRF are shown as function of τ -water of clouds, LWP and CFC changes. It is evident that the strong relationships in AMJ and JAS between total CRF, COT, and LWP are more important in modulating radiation than CFC. While SW CRF cooling dominates over LW CRF warming, CFC changes modulate mainly the LW portion of cloud radiation in both seasons. Arctic regionality emerges from the clustering of the regions, especially in AMJ and to a lesser extent in JAS. In the last two decades the net energy radiating at the surface due to clouds is decreasing.

350
Clouds cool the surface when the upwelling SW energy dominates the downwelling LW radiation. This is the case when clouds become optically denser and hence more reflective.
Those regions characterised by a darkening surface undergo an increase in SW reflection, leading to an increasing cooling by clouds (∆CRF<0). This takes place over the periennal sea ice zone (Beaufort, Laptev and East Siberian Seas), where a CRF decrease at a rate of −1-2 W m −2 is associated with greater cloudiness in AMJ and increasing τ -water in JAS. Quantitatively,

355
with values of ∆CRF Total = −1.4 W m −2 and ∆CF = 3.03 %, we obtain the total long-term sensitivity ∆CRF Total /∆CF = −0.48 W m −2 % −1 over the Beaufort Sea in AMJ. The sensitivities of the SW and LW parts of CRF amount to −0.56 and +0.84 W m −2 % −1 . Although averaged over one multi-year season only, our estimation is in line with measurements reported at the same location during the SHEBA campaign (Shupe and Intrieri, 2004). The SHEBA sensitivity of ∂CRF LW /∂CF = 0.65 W m −2 % −1 was seen to offset the SW for most of the year (with ∂CRF SW /∂CF ∈ [0, 1] W m −2 % −1 ), thereby warming the 360 surface while cloud cooling took place only in midsummer months. Accordingly, we report a net total (SW+LW) sensitivity of −0.13 W m −2 % −1 in JAS, meaning that the SW cooling takes over LW warming during the Arctic JAS in the record. The greenhouse effect from increased CFC in AMJ over these regions is not directly linked to the retreat of sea ice, the onset of which is in late May (Smith et al., 2020). Rather it is attributed to the enhanced convergence of atmospheric water content originating from open Arctic oceans during years with anomalously low sea ice extent (Kapsch et al., 2013). Our results imply 365 that this mechanism is not only evident in the year-to-year variability of exceptional sea ice lows, but is also a long-term component at decadal time scales, during which atmosphere-ocean coupling effects are predominant.
With the sole exception of the East Siberian Sea in JAS where τ -water of clouds grows in spite of a lower content of liquid water (∆r eff ≈ +0.3%, see Fig. C1), any positive τ -water trend corresponds to LWP changes for both seasons. Although not surprising, we note that the AMJ changes in CRF do not correlate with either LWP or COT. In the JAS months, however, 370 larger cloud optical densities and LWPs are matched by a decrease in CRF at the surface. This is the combined outcome of two effects: more insolation in the JAS months results in a more efficient SW scattering by cloud droplets and the darkening of the surface lowers the LWP threshold necessary for the CRF SW to dominate CRF LW . Excluding Barents Sea, the variance of ∆CRF during AMJ is narrower (−4.2 to +0.9 W m −2 ) than during JAS (−6 to +0.4 W m −2 ). This is evidence for the importance of radiance from the underlying surface, which is larger in AMJ than in JAS.

375
In addition to loss of clouds, the extended ice melt in Greenland is also known to be reinforced by low altitude liquid water clouds having sufficient opacity to allow the enhancement of downwelling LW flux but also being optically thin enough to allow a significant fraction of the SW flux to pass. This results in the surface being warmed (Bennartz et al., 2013). Such clouds occur within the range of LWP between 10 g m −2 and 60 g m −2 corresponding to the range of enhanced CRF at BOA. Fig. 8 shows that the increase of cloud τ -water and LWP over Greenland is among the largest throughout the whole Arctic for 380 both seasons. We attribute a reduction in melting to the increased wetting of clouds (∆LWP 40% at current typical LWP values). Similarly, during AMJ, Greenland and the Canadian Archipelago exhibit a rather small ∆CRF along with a large increase in τ -water and R TOA 560 . Enhanced surface cooling by clouds is thus less efficient. These regions are predominantly land or land/water mixtures, accordingly modulating the SW-to-LW flux balance.
Advances in observational techniques and process-level research are needed to assess unambiguously the relative roles of 385 temperature T and atmospheric particulate matter in determining cloud thermodynamic changes. In the absence of a systematic, pan-Arctic, aerosol indirect effect due to decreasing trends of ice or cloud condensation nuclei (IN/CCN), higher condensation rates (i.e. positive LWP trends) of small-sized cloud droplets can only nucleate and grow by a combination of changes in Arctic boundary layer depth within a saturated air volume. Different temperature regimes influence CA changing the τ − r eff −LWP relationship (Tselioudis et al., 1992) and favour droplet growth over condensation rates and vice versa (Lohmann et al., 2000).

390
To this end, the driver of, mostly decreasing, r eff trends (see bottom plot of Fig. C1) remains unclear. r eff size spectrum is modulated by the amount of water vapor and available particulate. While model and satellite data show a general moistening of the Arctic (Rinke et al., 2019;Boisvert and Stroeve, 2015), local on-ground (Graßl and Ritter, 2019) and pan-Arctic satellite (Linlu et al., 2021) evidence of a decrease in total aerosol burden is growing. However, IN/CCN can not be directly inferred from changes of column-integrated extinction of total aerosol load. Assuming a CCN decrease is in contradiction with the r eff 395 reduction via the Twomey effect. Alternatively, we speculate that the change in size spectrum or aerosol type might lead to optimal IN/CCN size and hygroscopicity, although the total aerosol amount has decreased. Satellite-derived single r eff values are only representative of the droplet/crystal population at a level of ≈ 1-τ from the cloud top (Platnick, 2000). We recommend that the available and relevant spectral observations are exploited King and Vaughan, 2012) to generate a pan-Arctic picture of in-cloud r eff (z) profiles, which would optimally complement surveys based on spaceborne 400 active techniques (Chan and Comiso, 2013;Matus and L'Ecuyer, 2017). r eff (z) profiles, together with aerosol speciation at high latitudes (Schmale et al., 2021) and cloud bases (Lelli and Vountas, 2018), are essential in two ways. First, they constrain IN/CCN activation, supersaturation and, therefore, cloud particle number concentrations (Zheng et al., 2015;Grosvenor et al., 2018). Second, cloud fields will be more accurately separated according to their phase (liquid, ice and mixed-phase) and layering (low, mid, high-level and multi-layered). We consider our results as upper bounds and more vertical resolution will 405 improve our understanding of the evolution of clouds in the Arctic.
From a modelling standpoint, we can validate past results (Morrison et al., 2019), for which the increase in cloud τ -water and LWP are projected to extend well beyond the middle of the present century. Constraining the cloud microphysics and thermodynamic phase will not only be crucial to project future Greenland melting (Hofer et al., 2019) but also to assess the sign and strengths of total cloud feedbacks (Gettelman and Sherwood, 2016). Given the actual and future Arctic temperatures, 410 ice will be increasingly depleted. Hence, τ -water and LWP will increasingly determine net cloud feedbacks (Bjordal et al., 2020). When the cloud ice phase turns to water a negative feedback is expected due to the offsetting of LW by SW. If climate models do not correctly capture this, i.e they do not incorporate more supercooled liquid than mixed-phase clouds (Lohmann, 2002), unrealistically large amounts of ice result, effectively reversing the sign of the net cloud feedback. We consider that this is one reason, which may explain in part the discrepancy between the atmospheric components (CAM) of the Community Earth that implies an increasing amount of supercooled liquid cloud droplets. The higher reflectance of clouds results in a more negative radiative forcing at the surface, thereby locally dampening Arctic Amplification. However, cooling by clouds implies the strengthening of the meridional temperature gradient. We expect this will lead to increase the inflow of warmer and moister 425 air masses from the lower latitudes into the Arctic climate. This may then either further decrease Arctic Amplification by generating more supercooled liquid water cloud, or possibly enhance Arctic Amplification by the increase input of warmer air, or some combination of the two. Future model projections of the Arctic climate must take into account these effects to accurately predict the impact of anthropogenic emissions of greenhouse gases and short-lived climate pollutants. Table A2 shows that overpass time, swath and footprint size differ among the sensors used in this work. These sensors are payloads on satellites which fly in sun synchronous orbits having different equator crossing times. Errors in the R TOA λ in the Arctic arising from the 30-minutes time lag are considered negligible for averaged R TOA λ . Monthly aggregation leads to higher means for finer spatially-resolved instruments than otherwise. Thus, intra-sensor radiometric R TOA λ harmonization is a prerequisite for the creation of calibrated time series and the detection of trends. Different application-dependent approaches have been 435 already employed. Krijger et al. (2007) derives gain correction factors based on the number of cloud-free scenes as function of spatial resolution for maximization of usable trace-gas retrievals. Tilstra et al. (2012) separates the influence of scattering geometry and cloud occurrence to correct SCIAMACHY reflectances for the computation of the aerosol absorbing index at UV wavelengths. Both approaches are not suited for our goal. The former aims at the removal of the influence of clouds, which are a primary component of the Arctic environment. The latter examines instrumental performance in a spectral region that 440 is not of direct interest as a result of potential radiometric degradation of sensors and of higher sensitivity to aerosols, whose radiative effects are comparatively small in the troposphere. Conversely, Hilboll et al. (2013) elaborate a method to explicitly take into account the difference in the ground pixel size and spatial misalignment across sensors. This is achieved by projecting the orbit of one instrument onto that of a second instrument. In our case, we select SCIAMACHY as reference sensor due to its well-calibrated spectral behaviour and because it overlaps with both GOME and GOME-2A. A conservative area-weighted 445 remapping scheme Jones (1999) is employed to derive the factor matrix transforming GOME-2A reflectances as they were measured by SCIAMACHY. Due to the frequent overlaps at high latitudes, only those GOME-2A orbits closest in time to SCIAMACHY are remapped. To extend the time series beyond the loss of Envisat in April 8th, 2012, full SCIAMACHY geolocations, comprising 431 orbits per month, have been used as target tessellation for the rest of the GOME-2A record. The downside of mimicking SCIAMACHY orbits, due to its design of alternating nadir and limb swath states, is the reduction of 450 the GOME-2A sampling rate. This is compensated for in part by the inherently different cross-swath viewing geometries and changes in illumination. GOME projection onto SCIAMACHY has not been implemented. Not only the two sensors overlap for a limited period of six months, but the relatively low sampling rate of GOME would have resulted in suboptimal statistics, even at a monthly scale. Validation has shown that GOME R TOA λ are consistent with those of SCIAMACHY (see Fig. 3 in main text). Remaining intra-sensor inconsistencies that cannot be compensated for, such as changes due to the dynamic radiometric 455 response over dark-to-bright surfaces, will be eventually accounted for by the trend model.
We tested the assumption that bidirectional surface effects do not introduce error in the detection of the temporal trends of R TOA λ by inspecting monthly distributions of the scattering angle throughout the record, separately for each sensor. This is needed because R TOA λ is, by definition, a directional quantity and depends on the scattering geometry. It has been found that the monthly data aggregation of all individual line-of-sights almost fully cover the hemispheric solid angle, incorporating The µ  Block bootstrap resampling (Efron and Tibshirani, 1993), belonging to the group of nonparametric methods, does not require prior knowledge of the analytical form of the underlying statistics of potentially non-normal data (Mudelsee, 2010). They rest on the block length of the effective independent random sample (Wilks, 1997, Eq. 19). An empirical sample distribution of the trend magnitude ω is then computed scrambling n times the blocks of the original record. The resulting empirical distribution 485 approximates the unknown ω probability density function. This allows to find the 2-σ ω interval needed for a confidence level at 95%. For all locations where the ratio |ω/σ ω | > 2, the trend magnitude ω exceeds natural variability and is termed statistically significant.
Appendix C: Uncertainty propagation in the cloud record and sensitivity The cloud data set is generated using an optimal estimation framework, which allows propagation of random and systematic 490 uncertainties into the pixel-based retrievals. Following Eqs. 2-5 in Stengel et al. (2017), for each location i at time t, we calculate the true variability σ true (i, t) and the uncertainty of the mean σ x (i, t) for the cloud property x from the mean of the squared pixel-based uncertainties σ 2 (i, t) and its standard deviation σ SD (i, t). Further, aggregation into monthly averages requires the uncertainty correlation c, or heterogeneity, relating σ SD (i, t) to σ true (i, t). Because c is not known beforehand, setting it to a fixed value is an arbitrary choice that does not account for the spatial and temporal relationship of algorithmic errors at pixel level throughout wide-scale cloud fields. Hence, we exploit the fact that σ SD → σ true when c → 1. This holds when the spatial sampling is the highest, thus we scale the number of successful retrievals of the cloud property x to c ∈ (0,1] and compute the c-dependent σ true (i, t) and σ x (i, t). Temporally, both σ true and σ x change as function of c. Seasonal trends of c reveal an overall increase of maximum 3% in AMJ and 1.9% in JAS over the Barents throughout the East Siberian Sea, whereas c over Greenland, Hudson Bay and the Canadian Archipelago exhibits a decrease of 0.6% in both seasons. This 500 translates into a change of ±0.5% and ±0.4% in σ true and σ x , respectively. With this approach, the clouds' heterogeneity of the monthly averages is related to retrieval errors predominantly in the spatial but not in the temporal dimension. Limited to an observational analysis of the cloud record, while uncritical for trend assessments only, σ x can be then successively used to label as meaningful those sensitivities of CRF to susceptible cloud property x, whose trend exceeds σ x .
Appendix D: Additional description of ozone trends 505 R TOA trends at 560 and 620 nm capture the Chappuis ozone absorption band having a broadband maximum centred about 602 nm and two wings stretching between 525 and 675 nm (Gorshelev et al., 2014). Analysing seasonal stratospheric and total column ozone, we are able to determine an effective modulation of R TOA trends by ozone. Ozone data in Fig 3 km is selected due to its highest sensitivity to stratospheric ozone concentrations, which peaks about that altitude. Ozone is produced in the tropics and circulation patterns transport it poleward. It is usually located above the tropopause and its concentrations are higher during winter months and lowest in summer months. Despite its high variability through the year, total ozone trends are generally small in the order of ±1%. Focusing on the Arctic, average total ozone is 353 DU and also 515 exhibits a distinct maximum in spring months and a minimum in summer months. The Arctic-wide trend of total ozone is positive by 3.9 DU (+1.1%) decade −1 , in line with global values.
Greater significant positive trends, ranging from +4 to +10% decade −1 , are found in stratospheric ozone. They are centred above Greenland and stretch out along the 75 • N parallel from the Greenland Sea through the Beaufort Sea in spring (AMJ) with a longer tongue over the Siberian Continent in summer (JAS). Contrasting the total with the stratospheric column yields 520 the influence of the tropospheric ozone only. For those locations where the trend in total ozone is absent but positive in the stratosphere, a negative tropospheric trend can be deduced. This mechanism is consistently found above 70 • N from the Canadian Archipelago through to the East Siberian Sea, irrespective of the season, together with the sustained positive trend above the Atlantic (the Greenland Sea), the neighbouring Barents Sea and the northern part of mainland Greenland (Gaudel et al., 2020). This reverses in a dipole fashion in JAS, when patterns of positive trends in total ozone are advected southward. In most eastern Arctic sectors (East Siberian, Laptev and Kara Seas) have a smaller contribution from ozone changes than the western sectors. This is consistent with a neutral ozone trend observed over these areas.
Finally, we speculate that a surface warming of the Arctic might inflate the tropopause, inducing the production of polar    Figure B1. Definition of the Arctic climate zones, identified by distinct geophysical settings, that will be used in this study to derive local trends of R TOA λ , cloud properties and forcing. The geographical subdivision follows that of Serreze and Barry (2014) and Wang and Key (2005a). at TOA and BOA.