Atmospheric trends over the Arctic Ocean in simulations from the Coordinated Regional Downscaling Experiment (CORDEX) and their driving GCMs

The Arctic Coordinated Regional Downscaling Experiment (Arctic-CORDEX) uses regional climate models (RCMs) to downscale selected Fifth Coupled Model Intercomparison Project simulations, allowing trend validation and projection on subregional scales. For 1986–2015, the CORDEX seasonal-average near-surface temperature (tas), wind speed (sfcWind), precipitation (pr) and snowfall (prsn) trends are generally consistent with analyses/observations for the Arctic Ocean regions considered. The projected Representative Concentration Pathway 8.5 (RCP8.5) 2016–2100 subregional annual tas trends range from 0.03 to 0.18 K/year. Projected annual pr and prsn trends have a large inter-model spread centered around approximately 5.0 × 10–8 mm/s/year and −5.0 × 10–8 mm/s/year, respectively, while projected sfcWind summer and winter trends range between 0.0 and 0.4 m/s/year. For all variables except prsn, and sometimes total precipitation, the driving general circulation model (GCM) dominates the trends, however there is a tendency for the GCMs to underestimate the sfcWind trends compared to the RCMs. Subtracting the Arctic-Ocean mean from subregional trends reveals a consistent, qualitative anomaly pattern in several variables and seasons characterized by greater-than or average trends in the central and Siberian Arctic Ocean and lesser or average trends in the Atlantic Sector and the Bering Sea, related to summer sea-ice trends. In particular, a strong proportional relationship exists between the summer sea-ice concentration and fall tas and sfcWind trend anomalies. The RCP4.5 annual, multi-model mean trends are 35–55% of the corresponding RCP8.5 trends for most variables and subregions.


Introduction
Significant near-surface temperature changes have occurred earlier in the Arctic than elsewhere on the globe and are expected to be more dramatic in the future (IPCC 2019(IPCC , 2021Post et al. 2019). Jansen et al. (2020) have compared modern trends with previous abrupt Arctic warmings and argue that climate models may even underestimate Arctic warming. Enhanced Arctic warming, referred to as "Arctic amplification" (e.g Manabe and Stouffer 1980;Koenigk et al. 2020) or more general "polar amplification" (e.g., Goosse et al. 2018) is linked to ice albedo and temperature feedbacks in addition to water vapor and cloud feedbacks (Cao et al. 2017), changes in atmospheric (Clark and Lee 2019) and oceanic (Beer et al. 2020;Koenigk and Brodeau 2014;Nummelin et al. 2017) heat transport and potential phytoplankton feedbacks (Park et al. 2015). Recent reviews and model analyses aimed at identifying the relative importance of Arctic amplification processes and feedbacks (Goose et al. 2018;Pithan and Mauritsen 2014;Stuecker et al. 2018;Zhang et al. 2018;Zhang et al. 2020) decompose the temperature feedback into separate contributions including a Planck feedback due to radiation changes caused by vertically uniform warming of the surface and troposphere and a lapse rate feedback due to vertically non-uniform warming.
These studies identified a relatively large and positive lapse rate feedback as the largest contributor to Arctic amplification, followed by a large positive surface albedo feedback and a smaller Planck feedback. However, the studies also highlight the difficulty of quantifying the feedback strength due to the non-linearity of many polar feedbacks and the dependence of feedback magnitude on location, season and climate state. This indicates that associated trends may vary significantly across regions. Satellite records of the recent warming in the Arctic also show a large and increasing spatial variation (Post et al. 2019), and Shu et al. (2020) suggest a link between the summer tendencies after 2000-including a faster decline in Arctic sea-ice extent-and spatiallyconstrained acceleration in the pace of melt which CMIP6 models had difficulty reproducing. Moreover, the Arctic climate appears to have entered an uncharacteristically steady anticyclonic regime since 1996 so that some climate indices no longer appear to be correlated as they were in the past (Proshutinsky et al. 2015).
Analysis of atmospheric trends over Arctic Ocean subregions is particularly relevant to ocean and ecosystem research. One of the key messages of recent Arctic ocean climate change and ecosystem assessments is the need for higher resolution basin-scale ocean-biogeochemistry models to provide locally applicable projections relevant for Arctic communities and management units (Steiner et al. 2015a;AMAP 2017AMAP , 2018a. These models require higher resolution atmospheric forcing such as is provided through the Coordinated Regional Downscaling Experiment (COR-DEX) project (Giorgi et al. 2009;Jones et al. 2011) to accurately simulate the impacts of climate change on the Arctic Ocean, sea ice and associated ecosystems. This is particularly important for coastal regions, where local communities access marine resources. The ocean and sea ice are highly responsive to the trends and patterns of key atmospheric forcing variables, air temperature, total precipitation, snowfall and wind speed through warming, stratification changes and radiation and momentum transfer, for example. These changes are expected to have important direct and indirect effects on the associated ecosystems. Hence, a large part of the uncertainty in the ocean and ecosystem model response is driven by uncertainties in atmospheric forcing. A thorough analysis of these uncertainties is therefore paramount in assessing the trends and projections of the Arctic marine environment and ecosystem.
The human population in the Arctic depends strongly on the natural environment for subsistence, and, along with the flora and fauna are very vulnerable to climate-related changes such as strong warming, sea-ice loss and extreme weather events (AMAP 2017(AMAP , 2018aIPCC 2019;ICC-Alaska 2015;ICC-Canada 2008). Recent large changes in temperature and sea ice have already been suggested to cause transformations in local ecosystems (e.g., Huntington et al. 2020) and shifts in the distribution of key Arctic marine species (e.g., Steiner et al. 2019) mammal feeding and migration habits (Harwood et al. 2015;Lotze et al. 2019). Projected changes into the future will likely further impact sea-ice (Tedesco et al. 2019;Lannuzel et al. 2020), pelagic (Kortch et al. 2012) and benthic (Renault et al. 2019) ecosystems and fisheries potential (Cheung et al. 2009;Tai et al. 2019). Changes in the Arctic climate may also strongly influence climate at lower latitudes (Koenigk et al. 2020 and references therein). Hence, a better understanding of the expected changes through historical simulations and climate projections is crucial.
Global climate models provide a general idea of largescale trends, however, their relatively coarse horizontal and vertical resolution restricts the ability to resolve atmospheric and ocean processes on a scale of regional relevance. Consequently, the reports listed above contain very little information on subregional differences in trends within the Arctic. The strong dependence of the polar climate on sea ice, which is challenging for coarse-resolution climate models, also adds considerable uncertainty to model projections in the region (e.g., Karlsson et al. 2013). In addition, poor observational data coverage in the Arctic, combined with large natural variability, makes validation of models challenging, particularly on a subregional scale. So, while Arctic climate change is undoubtedly complex, there is a pressing need for more detailed climate projections with realistic uncertainties to inform policy decisions and drive localized process models. As this change is not uniform across the Arctic, an evaluation of subregional trends allows assessment of potentially different impacts and/or response time scales for different communities around the Arctic. It is also important to understand the degree to which the subregional trend patterns simulated in the driving global models are retained in the regional downscaling. Such consistency would allow us to assess the uncertainty in trend projections from global models while using individual downscaling models to force basin-scale ocean ecosystem models and assess local changes of importance to communities.
The Arctic-CORDEX project includes nine atmospheric regional climate model (RCM) configurations that were run for the historical period and the Intergovernmental Panel on Climate Change (IPCC) RCP8.5 scenario [littleto-no action to curb greenhouse gas (GHG) emissions, IPCC (2014)] on high-resolution Arctic domains (ARC-22 and ARC-44, with resolution 0.22 and 0.44 degrees, respectively). Six of these were also run for an intermediate scenario, RCP4.5 and an additional configuration has data available for RCP8.5 only. These simulations, summarized in Table 1, are particularly well-suited to a detailed investigation of Arctic climate change and uncertainty because they include cases of a single RCM driven by several global general circulation models (GCMs), as well as different RCMs driven by the same GCM, allowing the effects of the driving simulation and downscaling model to be compared.
Daily precipitation and temperature indices from the CORDEX historical Arctic simulations for 1980-2004 have been assessed over Canadian land by Diaconescu et al. (2018) who found that the RCMs generally simulated mean temperatures and warm extremes well, but were less successful at simulating indices related to cold extremes and precipitation. They also found that the ARC-22 version of CanRCM4 performed somewhat better over land than the lower-resolution configuration. The four RCA4 configurations have been analyzed by Koenigk et al. (2015) who found that while, unsurprisingly, the RCM simulations were strongly influenced by the driving GCM, the projected tas changes were significantly different over most of the Arctic Ocean and there were some interesting differences such as a reduced winter 2 m temperature bias over the Arctic Ocean for the RCM ensemble mean compared to that of the GCMs (with respect to ERA-Interim) and an opposite future cloud cover response under RCP8.5. They also reported that the Arctic average precipitation and temperature response under RCP8.5 and RCP4.5 are linearly related. Future projections of Arctic cyclone frequency and properties were investigated using a similar set of simulations by Akperov et al (2019). They found that, while the lateral boundary conditions from the driving GCM had a stronger influence on the projections than the RCMs, the RCM simulations showed opposite changes in cyclone frequency to the driving GCMs over the central Arctic Ocean and East Siberian Sea in winter and most of the central Arctic in summer.
Linear trends of key variables provide a compact representation of the lowest-order subregional climate change and variability on a multi-decadal time scale. This is useful for validation of the performance of this suite of simulations over the recent past, as well as to summarize future projections, though sensitivity to choice of endpoints must be considered. While the linear approximation becomes less valid with increasing timescale, the confounding influence of natural internal variability decreases, allowing the forced component of the climate-change pattern to dominate. Here, the recent historical trends for 2 m air temperature (tas), total precipitation (pr), snowfall (prsn) and 10 m wind speed (sfcWind) over Arctic ocean basin subregions for the COR-DEX models and their driving GCM simulations are validated with respect to selected analyses and observations. The relative influence of downscaling and boundary conditions on these trends is also explored. Future projected trends and uncertainty for these subregions, according to scenarios RCP8.5 and RCP4.5, is presented, as well as their relationship to trends in projected sea-ice concentration, and possible implications for Arctic Ocean ecosystems are discussed.

Models and data
The Arctic CORDEX simulations use various atmospheric RCMs with lateral and surface (ocean, sea-ice and land) boundary conditions provided by selected GCMs. The historical, RCP4.5 and RCP8.5 scenario CanRCM4 simulations  (Brodeau et al. 2009;Dussin et al. 2016) and GISTEMP v4 (GISTEMP Team 2020; Lenssen et al. 2019) data were included and for pr, the GPCP data set (Adler et al. 2003) was also used. The reanalyses and GISTEMP data are complementary, the former using model assimilation to produce a high-resolution (especially in the case of ERA5) global data set from sparse atmospheric data, and the latter estimating surface air temperatures from relatively well-measured sea-surface temperatures and combining them with nearsurface air temperatures from land stations using binning by equal-area squares and sophisticated quality control procedures (Hansen et al. 2010). ERA5's much higher resolution makes it more suitable for validation of RCMs, but GISTEMP is more directly tied to the observations. Since this work is partially motivated by the desire to optimize the transition between historical and future atmospheric forcing of ocean/ice models, temperature trends from the ERA-Interim (Berrisford et al. 2011) based ocean-forcing data set, DFS, which includes bias-corrections in the Arctic to accommodate more realistic sea ice, were also included.
Performance of ERA5, the highest-resolution data set available, and thus most suitable for comparison with COR-DEX simulations, is mixed over the Arctic Ocean; it outperformed the ERA-Interim, JRA-55, CFSv2 and MERRA-2 analyses when compared to temperature, wind speed and specific humidity buoy data from the Fram Strait, late-summer 2017 (Graham et al. 2019b) and had the best spring and summer downward radiative fluxes over sea ice but worst temperature biases during the Norwegian Young Sea Ice Campaign (N-ICE2015) (Graham et al. 2019a); in Wang et al. (2019) ERA5 was seen to correct the known Arctic snowfall/total-precipitation ratio bias of ERA-Interim but showed larger temperature biases than ERA-Interim compared to Arctic sea ice buoy data for 2010-2016. The inclusion of DFS and GISTEMP, as well as the other analyses, in our tas validation provides some confidence in our conclusions regarding the trends. Likewise, while calculation of trends from observational data should be treated with caution in data-sparse areas, the inclusion of GPCP in our precipitation validation as well as other reanalyses for all variables gives an indication of the robustness of these results.

Analysis
The subregions considered in this report, shown in Fig. 1, are based on the Arctic Ocean basins 1-14 (Carmack and Wassmann 2006;Matrai et al. 2013). The masks for these regions, as well as the ARC-44 CORDEX simulations, CanESM2 data and analyses/observations were all remapped onto the CORDEX ARC-22 grid. For the masks, nearest neighbor interpolation was used, followed by manual editing for consistency with the CORDEX ARC-22 land-sea mask (Online Resource 1). For the model and observation-based data, bilinear interpolation was used.
Seasonal and annual linear tas, pr, prsn and sfcWind trends for the historical period, 1986-2015, were calculated for each subregion from the RCM simulations and compared with corresponding trends from the analyses and observations, as well as available GCM drivers. Estimates of the effect of natural variability on the historical tas, pr and sfcWind annual trends, as well as RCP85 tas trends, were obtained as 95% confidence intervals, computed as twice the standard deviation of the corresponding trends from the 50-member large ensemble of the CanESM2 GCM, in addition to those computed from the CORDEX suite. Spatial tas trend patterns over the Arctic Ocean were also compared with ERA5 using Taylor diagrams.
Projected 2016-2100 subregional linear trends were then calculated for the CORDEX and driving models with available RCP8.5 simulations, including HIRHAM5-v2, ESM-LR which missed the validation step, and RCP4.5 simulations. Trend anomalies for the subregions were defined for each simulation as the residual when its Arctic-Ocean mean trend is subtracted, i.e., for each model, m , and subregion, s , the trend anomaly, A ms , is given by A ms = T ms − P m , where T ms is the raw subregional trend and P m is the pan-Arctic-Ocean trend, defined as the trend of the variable averaged over all subregions in Fig. 1. Confidence intervals for the corresponding driving-model trends were calculated based on the residual variance of the single ensemble members, since their purpose is to evaluate the statistical significance of trend differences between the driving models and those of the downscaling RCMs for a given realization.
Two types of multi-model mean trends were calculated, representing the two extreme sets of assumptions regarding the importance of the driving GCM in determining the trends. The straight multi-model average weights the trend of each model equally, based on the unlikely assumption that models with the same driving GCM contribute equally to the trend variability. The "weighted mean" weights each model inversely by how many models share its driving GCM according to the extreme assumption that the driving model is entirely responsible for the trend variation. The best possible estimate of the mean available from these simulations should lie between these values.
For insight into the geographic variation of the long-term trends, least-squares fits to simple linear models of the trend patterns were also performed.

Results
The annual average recent simulated historical and RCP85 projected regional linear trends from the CORDEX suite for tas, pr, prsn and sfcWind are presented below. Selected seasons are also shown for some variables, along with subregional trend anomaly patterns and a comparison between RCP8.5 and RCP4.5 annual average trends as well as some corresponding trends for the driving GCMs. The mean and uncertainty of all of the seasonal tas, pr and sfcWind seasonal linear trends for 1986-2015, as well as RCP8.5 and RCP4.5 future projections are given in Table S1 of the supplementary material.

Historical trends
3.1.1 2 m air temperature Figure 2 shows the regional linear trends in annual average tas for 1986-2015 for the Arctic CORDEX historical (1986-2005) + RCP8.5 (2006-2015 simulations. The trends from the GCM ensemble members used to drive the COR-DEX models (except EC-Earth_r3, which was not available on the Earth System Grid Federation (ESGF) server) and the observation-based data sets are also shown. Note that the RCPs are very similar in the early part of this century, so these results should not depend strongly on the choice of the RCP8.5 scenario for the last decade. Also shown are the mean of the CORDEX models and 95% confidence interval error bars estimated (as double the sample standard deviation, 2σ) from the CanESM2 50-member large ensemble in black, and the CORDEX suite, in grey. The two error bars are usually similar, suggesting that natural variability likely dominates the inter-model spread in most subregions, except possibly 5, 12 and 13.
The CORDEX models show 1986-2015 annual tas trends ranging from about 0.0 to 0.3 K/year over the Arctic Ocean subregions. The straight model means are slightly larger than the weighted means, in most cases, but both are between approximately 0.05 and 0.15 K/year. In most subregions the observation-based trends are not significantly different from the CORDEX mean at the 5% level using either set of error bars. The "observational trends" are generally within the 95% confidence intervals of the simulated trends and there are no cases where all observational estimates lie outside them. This suggests that the trends from the suite of simulations are reasonably consistent with those of the observation-based data sets for the recent past. It is unsurprising that the magnitudes of the trends from models driven by the same GCM ensemble member (indicated with the same color) tend to be quite close, while the simulations using the same RCM driven by different GCMs (indicated with the same shape) tend to differ more. As a further illustration of the dominance of the driving models, the CORDEX trends are quite close to the corresponding driving-model trends. For example, CanRCM4 (which is spectrally nudged), and to a lesser extent, RCA4 (with no nudging) subregional trends are very close to the corresponding values from the CanESM2 driving ensemble member. There is also very little effect of downscaling on the pan-Arctic Ocean mean trend. The RCA4_NORESM1 and HIRHAM5-v1_EC-EARTH-r3 simulations are generally low outliers while the simulations driven by CanESM2 and EC-EARTH-r12 1 3 tend to be on the higher side. The subregional variation is considerably less than the variation among models for this time period and there does not appear to be a significant qualitative pattern. All observation-based data sets (except MERRA-2 in the Canadian Archipelago) and all but one of the models (RCA4_NorESM1 in the Southern Beaufort Sea) give a positive trend for all regions. Figure 3 shows the corresponding seasonal trends for the CORDEX simulations and observation-based data. Winter (DJF) was computed using the December of the previous year. Summer (JJA) generally has the weakest trends for both models and observation-based data. Autumn and winter have generally strong trends and large inter-model spread, with several models exhibiting negative trends in Baffin Bay (region 3) and one model, RCA4_NorESM1, exhibiting negative trends in all three subregions covering the Canadian Arctic Archipelago (3,4 and 5) in winter. In all seasons the "observed" trends are within, or close to, the spread of the CORDEX trends, in most cases. Note that since the free historical simulations were used, rather than those forced by reanalysis data, as were analyzed in Glisan et al. (2019) and Akperov et al. (2018) they are not expected to reproduce the "realization" represented by the observations.
The spatial correlation of the tas trends over the Arctic Ocean with ERA5 is between 0.3 and 0.7 for all models in both DJF and JJA (Online Resource 2, Figure S1). Interestingly, the RCA4 and RCA4-SN simulations driven by EC-EARTH_r12 display an approximately 50% larger standard deviation of the trend pattern than ERA5, while HIRHAM5-v1, driven by EC-EARTH_r3 displays close to a 50% smaller standard deviation in both seasons. Care must be taken not to over-interpret the correlations as they are much more susceptible to natural variability than those of the climatologies, which have much larger correlations with ERA5 (Online Resource 2, Figure S2). For example, the correlation between the ERA5 annual trend pattern and that of the very similar DFS data is only 0.7. As another illustration of this sensitivity, Figure S3 in Online Resource 2 shows winter and summer Taylor diagrams of ERA-5 with time-shifted versions of itself. It takes only a shift of eight to ten years, which should retain most of the forced trend pattern but change the representation of the natural variability, for the spatial correlations to be reduced to approximately 0.5. The main message to take from the correlations in these plots is that all models have spatial correlations of the trend patterns with ERA5 over the Arctic Ocean consistent with having a similar forced trend pattern, as they correlate as well with ERA5 as approximately 10-year shifted versions of ERA5.

Precipitation
The 1986-2015 historical + RCP8.5 subregional total precipitation (rain and snow) annual linear trends are shown in Fig. 4. The multi-model means are small and positive for all subregions with both methods, straight and weighted (not shown), of calculating the multi-model mean giving similar results. The trends are not likely significantly different Fig. 2 The subregional (see Fig. 1) 1986-2015 annual average tas trends (K/y) for the CORDEX historical (1986( -2005( ) + RCP8.5 (2006( -2015 simulations (open shapes) and comparison data sets (asterisks) with the mean of the CORDEX models (large filled circles), the driving GCM ensemble members (small filled circles) and 2σ error bars estimated from the CanESM2 historical large ensemble (black) and CORDEX suite (grey) Fig. 3 The subregional (see Fig. 1) seasonal 1986-2015 tas trends (K/y) for the CORDEX historical + RCP8.5 simulations, the mean and weighted mean (described in Sect. 2.2) of the CORDEX models (filled circles), 2σ error bars estimated from the CORDEX suite and comparison data sets with the DJF time series running from December 1985 to February 2015

Fig. 4
The subregional (see Fig. 1) 1986-2015 annual average pr trends (mm/s/y) for the CORDEX historical (1986( -2005( ) + RCP8.5 (2006( -2015 simulations (open shapes) and observations/analyses (asterisks) with the mean of the CORDEX models (large filled circles), the driving CanESM2 ensemble member (small filled circles) and 2σ error bars estimated from the CanESM2 historical large ensemble (black) and CORDEX suite (grey) from zero for most subregions, based on the CORDEX and CanESM2 large-ensemble error bars. As with temperature, the CORDEX error bars are generally comparable to the natural variability estimated from the CanESM2 large ensemble and the majority of observations and analyses are within both sets of error bars, in most cases, except for the northern Beaufort Sea. There is a general tendency for the simulated trends to be larger than the "observed" trends. The driving GCM appears to be a major factor in the trends, with some exceptions. In particular, there is a relatively large spread in the ESM-LR-driven model pr trends in subregions 2 and 3, surrounding Greenland. The downscaling RCMs have very little effect on the pan-Arctic-Ocean mean trends.
The spread in simulated pr trends is generally greater in summer, where the ERA5 trend is consistently lower than the majority of the models ( Figure S4). In the other seasons, the ERA5 pr trends generally fall within the model 95% confidence intervals, which encompasses zero for most subregions and seasons. Note that, since precipitation rates in the Arctic are very low, even small trends can have important impacts.
The 1986-2015 snowfall trends ( Fig. 5) are small, with the mean trends generally negative or close to zero and surprisingly, given the complicated interplay of variables involved in determining snowfall, close to the "observed" trends, especially ERA5, in most cases. The negative trends are generally due to an increase in the pr fraction falling as rain, rather than a decrease in precipitation, since total precipitation trends are generally positive (Fig. 4). In all subareas, the ERA5 trends are within the model spread error bars, while JRA-55 and MERRA-2 are sometimes below. In the case of MERRA-2, there is a general positive bias with respect to GPCP in total precipitation over highlatitude oceans (Gelaro et al. 2017) which would translate into a stronger negative trend in prsn with warming. The error bars encompass zero for all subregions except the Bering Sea, the most southerly subregion considered. The CORDEX trends are often close to those of their driving GCM, especially for those driven by CanESM2. There are some exceptions, such as EC-Earth_r12-driven trends in subregions 2, 13 and 14 and ESM-LR-driven RRCM in subregion 2. The possible significance of downscaling effects on snowfall trends in these regions is investigated further in Sect. 3.2.2 using the future RCP8.5 projections, which have a stronger "signal to noise" ratio.
Zhang et al (2019) find that over the Arctic ocean, reductions in snow cover fraction over sea ice, and sea ice extent appear to contribute equally to the Arctic albedo decline with each of similar magnitude to reductions in terrestrial snow cover north of 60 N. The decrease in snow cover fraction is attributed primarily to the increase in surface air temperature, followed by declining snowfall. Consistent with our pan-Arctic reanalysis and most model results, Zhang et al (2019) indicate that, although total precipitation has increased, Arctic snowfall is reduced in all of their analyzed data sets due to warming. The 1986-2015 annual average prsn subregional trends in mm/s/y water equivalent, symbols as in Fig. 4

Winds
The historical simulated annual wind-speed trends, shown in Fig. 6, are generally small and span zero for all subregions except 12 and 13, where they are all positive and, in the case of simulations driven by EC-EARTH, quite large (up to 0.05 m/s/year). ERA5 and JRA-55 are also positive for those subregions, but Merra-2 displays slightly negative trends there. As with the other variables, the trends from the analyses are within the 95% confidence intervals around the multi-model means based on the CanESM2 large ensemble for most subregions, with MERRA-2 and JRA-55 as large negative outliers in subregions 4 and 5,10, respectively. Subregions 4 and 5 are likely particularly affected by the lower resolution of these data sets, as they lack large areas of open ocean. The 2σ error bars from the CORDEX suite are similar to those from CanESM2 in most cases. Exceptions are the East Siberian Sea, with CORDEX standard deviation about half that of CanESM, and the Kara and Barents where the EC-EARTH-r12 driven model trends (purple) are high outliers, contributing to a much larger standard deviation for the CORDEX simulations. The straight and weighted multi-model means are very close to each other and positive or near zero for all subregions. While the sfcWind trends tend to group quite closely by driving model, the actual driving model trends are often much farther from the corresponding RCM trends, suggesting a systematic effect of downscaling alone, even over ocean areas. This is further investigated using the RCP8.5 simulations in Sect. 3.2.3. Table 2a shows the tas, pr and sfcWind winter and summer and prsn annual subregional trends for the RCP8.5 "highemission" scenario Arctic CORDEX simulations. Particular features are discussed below. Since all simulations were close-to or within the estimated natural variability of the ERA5 trends and there were no consistently strong outliers in the Taylor diagrams, all simulations were included in the calculations for the future projections. An additional simulation available for RCP8.5, HIRHAM5-v2_ESM-LR, was also included. While that configuration could not be assessed in the validation step, its components (a different version in the case of the RCM, HIRHAM5-v1_EC-EARTH) were both represented there in different configurations. the CanESM2 large ensemble, indicating the likelihood of systematic trend differences between models. The colored lines give the means of trends in tas averaged over all the ocean subregions. For subregions 2, 13 and 14, representing the Atlantic in and outflow, and subregion 7, the Bering Sea (Pacific inflow) the subregional trends in all simulations are less than the mean trend of "Arctic-Ocean average" tas, a fairly robust indication that, in spite of the large uncertainty in the absolute trends, these areas experience weaker than average tas trends under this scenario. Similarly, subregion 1 (central Arctic) and subregions 6, 8, 9, 11 and 12 (Canadian Basin, Chukchi and Kara Seas) have consistently larger than average tas trends under RCP8.5. Taking anomalies of the subregional trends with respect to their Arctic Ocean averages to remove the influence of the pan-Arctic mean on the variation, as shown in the bottom plot of Fig. 7, reveals the consistency of the spatial trend pattern mentioned above more clearly. It also highlights the relatively small effect of the downscaling RCMs on the tas trend pattern at this scale. In most cases, the inter-model spread of the anomalies taken this way is a fraction of the spread of the raw trends, confirming that the pan-Arctic difference in the particular driving GCM simulation is a major contributor to the spread. For subregion 7, the Bering Sea, the trend anomaly spread is double that of the raw trend, a notable exception to this. The seasonal RCP8.5 tas trends in Fig. 8 indicate that the future DJF trends tend to be at least 3 × higher than the corresponding JJA trends with a much larger spread, similar to the historical trends. The subregional anomaly pattern in SON and DJF indicates generally smaller trends in the Atlantic sector, where there is less sea-ice cover, as well as the more southerly Bering Sea (region 7) and larger than average trends in the central and Siberian Arctic Ocean, similar to the annual average trends. In summer, the raw trends appear fairly uniform. The anomaly pattern, however, has several subregions where all models give trends above (12 and 13) or below (1, 6 and 9) the Arctic Ocean average. The spring and fall trends are generally intermediate in magnitude between those of summer and winter, but while the former has no discernible subregional trend anomaly pattern, the latter has, arguably, the clearest pattern of all. It is also interesting that, for SON, the overall pattern is amplified more in some models than others. In fact, the high degree to which the raw trends can be expressed as the sum of the pan-Arctic-Ocean trend and a scaled version of the weighted mean subregional anomaly pattern, A s , is given by the small symbols in Fig. 8d, which represent the residuals,R ms = A ms − f m A s , for model m and subregion s . f m is the scaling factor determined by a least squares fit of A s to the model's trend anomaly pattern,A ms . The largest residual magnitude is less than 33% of the corresponding raw trend and the vast majority are less than 20% of the corresponding raw trends. The tendency, noted in the annual average plots, for the spread in the anomalies to be smaller than in the raw trends generally holds for the individual seasons as well. As with the annual average, the Bering Sea (7) is an exception to this in all seasons. In particular, in fall the models agree amazingly well on the absolute trend in that subregion but disagree strongly on how much weaker it is than the Arctic-Ocean average. The winter anomaly trend spread for Baffin Bay is also somewhat larger than that of the absolute trends. The tendency, noted in the annual average plots, for the spread in the anomalies to be smaller than in the raw trends generally holds for the individual seasons as well. Figure 9a shows the annual linear total precipitation trends under RCP8.5 for the Arctic CORDEX simulations and available driving ensemble members, indicating a general, small increase across all subregions and seasons in this scenario except for a few instances in the RCA4, RCA4-SN and HIRHAM5-v2 simulations forced by ESM-LR. There is considerable inter-model spread in the projected total precipitation trends. In particular, the RRCM_ESM-LR simulation has very different behavior than all of the other simulations and ERA5. It is often an outlier in total precipitation with, most notably, an approximately five times stronger trend in Baffin Bay (3). It is interesting to note that the RCM simulations forced by the same ESM-LR ensemble member give quite different precipitation trends for some subregions, even more marked for the seasonal trends. There is, however, a Fig. 8 The 2016-2100 RCP8.5 seasonal linear tas trends and subregional trend anomalies (see Sect. 2.2) and residuals left when the Arctic-Ocean mean trends and least-squares fit of the weighted mean anomaly trend pattern are subtracted from the corresponding raw trends for SON (small symbols) general tendency for RCM total precipitation trends to be close to those of their driving models. The relatively large difference between EC-Earth_r12 pr annual historical trends and those of the corresponding RCMs are not reflected in the annual-average RCP8.5 case, though there are relatively large differences in the seasonal trends (not shown) for a number of subregions. Also, the JJA RCA4 precipitation trends (not shown) are often considerably higher than those of CanRCM4, forced by the same CanESM2 ensemble member and CanESM2, itself. This suggests, that downscaling models could have a relatively large effect on precipitation trends at a local and seasonal level, even though the annual CORDEX precipitation trends are generally within the 95% confidence interval of the corresponding GCM trends.

Precipitation
Considering the subregional trend anomalies, calculated with respect to the respective Arctic-Ocean average trends as was done for tas, we note that, in contrast to tas, the spread is similar to that of the absolute trends for many areas. This suggests that the overall variation in the pan-Arctic trend influences the spread very little in those subregions. There are also no Arctic Ocean subregions with all anomalies larger or smaller than zero, so there is no consensus on subregions likely to experience greater or less than average annual pr trends as there was for tas.
The snowfall trends in Fig. 9b illustrate the competing effects of increasing temperature (leading to snow-rain conversion) and increasing total precipitation. As a result, the decadal-scale behavior of the timeseries are not as linear for prsn as they are for other variables. In particular, HIRHAM5-v1_EC-EARTH shows an initial increase, related to increased total precipitation, followed by a decrease, as rising temperatures cause less precipitation to fall as snow, for many subregions. Unsurprisingly, the models tend to show decreases in snowfall for the subregions with larger fractions of their area farther south, as snow turns to rain, and trends close to zero for the more northerly subregions. Note, however, that the subregions with strong negative snowfall trends, the Greenland shelf (2) and Bering Sea (7) are not those with larger-than-average temperature trends. The multi-model mean trends are negative or near zero for all subregions. The CORDEX prsn trends are often significantly different from their driving-model trends. They are associated rather more by their downscaling RCM which appears to have a systematic influence. For example, the HIRHAM5-based simulations are near the top edge of the spread for all ocean subregions except moderate trends in the Bering Sea and east coast of Greenland and strong negative trends in 13 and 14 off Scandinavia. Looking at the subregional trend anomalies, in addition to the HIRHAM5-v2 and HIRHAM5-v1 simulations generally having similar values, the simulations using RCA4 and RCA4-SN tend to be clustered.

Winds
The winter and summer RCP8.5 sfcWind trends in Fig. 10 (top) are generally small and positive, with the exception of some, mostly ESM-LR-driven, simulations, which have very small negative trends in subregions 2, 3, 7, 13 and/or 14. The CORDEX trends are usually closely associated according to driving model, however, the driving model trends are often significantly lower than the corresponding CORDEX trends, further evidence of a systematic qualitative effect of downscaling alone. The annual monthly mean daily maximum wind speed (sfcWindmax) mean trends (not shown) are very similar to those of sfcWind, suggesting that brief high-wind events are not strongly preferentially affected by the simulated climate change. Unsurprisingly, the sfcWindmax trends have a greater spread.
The similarity between the RCP8.5 sfcWind trend patterns and some seasonal tas patterns is particularly apparent in the subregional trend anomaly plots (Figs. 10 and 8, bottom). In particular, the similarity between the SON tas trend anomaly pattern and that of sfcWind in JJA is uncanny. It is, at first glance, surprising that, in contrast to tas, for sfcWind this characteristic pattern is more evident in the JJA trends than in the DJF trends. The inter-model spread in the trend anomalies is generally less in summer than in winter, but the difference is less than that of tas. Also, as with tas, the inter-model spread in the anomalies is generally similar-to or smaller-than that of the absolute trends except for the Bering Sea (7) and Baffin Bay (3) in the winter.

Comparison of the RCP4.5 and RCP8.5 scenario simulations
Figures 11a-c show the subregional trends for 2016-2100 of the RCP4.5 simulations for tas, pr and sfcWind, respectively, as well as the corresponding RCP8.5 multi-model means for comparison. Note that the RCP4.5 means are based on a smaller number of simulations. For all three variables, the qualitative trend pattern is similar to, with about half the magnitude of the RCP8.5 pattern. The large degree to which the subregional trends scale linearly between scenarios is highlighted in Fig. 11d which shows the ratio of RCP4.5 multi-model means to the means of the RCP8.5 trends taken over the same subset of models (i.e., only model configurations that include both RCP4.5 and RCP 8.5 simulations are included) for tas, pr and sfcWind. For tas, the individual trend ratios for individual models are also plotted. These ratios are all between approximately 0.35 and 0.55 for tas and pr multi-model means. In fact, even the individual model ratios for tas only range from about 0.25 to 0.6 and the relative ratios of the different models are quite consistent, e.g., RCA4, ESM-LR has the smallest RCP4.5/RCP8.5 Fig. 10 The 2016-2100 RCP85 linear sfcWind a winter (DJF) and b summer (JJA) trends and subregional trend anomalies with filled symbols in the anomaly plots indicating negative absolute trends tas ratio for most subregions and RCA-4, NorESM1 and HIRHAM5-EC-EARTH r3 are usually in the high end of the range. Note, the relatively high ratios for subregion 13, where the strong changes seen in the historical period persist farther into the future projections. The tas and pr multimodel mean ratios are uncannily similar, differing by less than 0.1, usually much less, indicating that the linear relationship between total Arctic pr and tas change in the RCA4 models shown by Koenigk et al. (2015) also holds for the Arctic Ocean subregions, individually. The sfcWind ratios are also similar to those of tas and pr in subregions where the trends are not small compared to the spread. It would be interesting to see the degree to which this scaling remains in fully-coupled high-resolution models.

Historical trends and validation
For tas, the historical subregional trends are largely retained in the downscaling process, even for the simulations without spectral nudging. For pr, prsn and sfcWind, there are cases where downscaling introduces differences comparable to the differences between driving model trends, though for pr and sfcWind, the downscaled trends still tend to be closely grouped by driving models, suggesting that it is the downscaling itself, rather than the details of the RCM physics that is responsible, in most cases. Determining the consistency of the suite of free-running historical CORDEX simulations with observation-based data is an important step in evaluating the future trends projected by the models. Unfortunately, this cannot be done rigorously without sizable ensembles with which to estimate the natural variability of each model as, even for a perfect model, the observations would only be expected to represent one possible realization. In addition, the observation-based data currently provide only weak constraints on the historical trends because of the relatively short period of widespread Arctic measurements and high natural variability. To give some indication of the overall consistency between the subregional trends from this group of simulations and the recent historical observed trends, the natural variability was estimated from the CanESM2 historical large ensemble. It is not unlikely that this is a conservative estimate of the variability, as the necessary resolution limits and parametrizations of a GCM reduce the complexity of the system. As the historical CORDEX and observation-based trends are Fig. 11 The 2016-2100 RCP45 a tas, b pr and c sfcWind annual linear trends. d RCP4.5/RCP8.5 ratios of mean tas (red), pr, and sfcWind trends for model configurations with both scenarios available, as well as the individual model ratios for tas generally within or close to the estimated "error bars" of the straight and weighted multi-model means for all variables, the simulations and observations appear roughly consistent. It is also important to keep in mind that the analyses have relatively little observational constraint in the Arctic. In fact, Diaconescu et al. (2017) found that, compared to tas and pr indices computed from land station data, CanRCM4 outperformed ERA-Interim, though the lower resolution of ERA-Interim likely degraded its performance relative to station data more than would be expected of ERA5. Thus, the degree to which outliers may be less reliable cannot be determined with the analysis presented here, so all models were used for the future trend calculations.

Projected subregional trend patterns
The future subregional trends tend to group together much more by forcing GCM than downscaling RCM for all variables except snowfall and, in some cases, total precipitation. This grouping indicates that even though there can be statistically significant differences in climate-change signals between RCM simulations and their driving GCMs over the Arctic Ocean (Koenigk et al. 2015), the subregional trends are largely determined by the driving simulation. The grouping further suggests that trend uncertainties determined from multi-model and multi-ensemble global modelling efforts might translate to the downscaling model, at least for temperature, confirming the importance of considering the global-model natural variability when using atmospheric model data to force regional ocean or ecosystem models. However, it also suggests that systematic differences in downscaling RCM snowfall parametrizations may be a major source of uncertainty in that variable. The fact that the RCP8.5 tas trend standard deviation is generally much larger than the estimate of the 95% confidence limit of the estimated natural variability suggests that the systematic differences between model configurations dominate the intermodel spread in that case.
In the case of wind speed, even though the driving models appear to be the dominant influence, as demonstrated by the close grouping of CORDEX subregional trends by driving simulation, the sfcWind trends for the GCM simulations are generally lower than those of the corresponding CORDEX simulations. This suggests that smaller-scale wind features tend to increase or decrease along with the subregional-scale mean wind. It also indicates that mean wind speed-changes over the Arctic Ocean under climate change are likely underestimated by low-resolution general circulation models.
For total precipitation, the larger spread in the trends for different CORDEX simulations with the same driving ensemble member in some areas and seasons indicates, unsurprisingly, that precipitation trends are less tied to those of the driving model than those of temperature and wind speed. The larger RCP8.5 projected JJA precipitation increase in the RCA downscalings compared to the GCM simulations for the Arctic, in Koenigk et al. (2015) is generally reflected in our seasonal precipitation trends (not shown). CanRCM4 and RRCM tend to produce JJA precipitation trends closer to those of CanESM2 and ESM-LR, respectively, than the other CORDEX simulations, major exceptions to this being extremely large trends in RRCM around Greenland, which may be a topographical effect in the downscaling.
In spite of the large inter-model spread in individual trends, the consistency of the qualitative future Arctic subregional trend anomaly patterns between models, scenarios and, in some cases, even across different variables, is remarkable. This suggests that strong underlying mechanisms, common to all models, determine the trend anomaly pattern. This is illustrated in Table 2b, which summarizes the multi-model mean subregional trend anomalies for the RCP8.5 simulations. Cases where all simulations give values larger or smaller than their respective Arctic-Ocean means are indicated by dark orange or blue cells, respectively, and cases where at least 8 of 10 simulations give positive/negative local trend anomalies are indicated by light orange/blue. Note that for all variables except prsn, the Arctic-Ocean mean trends are positive, so negative trend anomalies usually mean that the subregion experiences weaker-than-average positive trends. A similar trend anomaly pattern applies to all except JJA tas (and possibly JJA pr, which is trivially consistent with no discernible pattern). It is characterized by greater-than or approximately average trends in the central Arctic subregion (1), lesser or average trends in much of the Atlantic Sector (2, 13 and 14), greater or average trends in the Siberian Arctic Ocean (8,9,10,11,12) and lesser or average trends in the Bering Sea (7). For RCP8.5 tas and sfcWind, subtracting the pan-Arctic trends reduces the intermodel spread considerably for most subregions, indicating that much of the subregional trend disagreement between models comes from their simulation of the overall Arctic response to the RCP8.5 forcing.
Recent analysis of the contributing feedback processes to the Arctic amplification signal indicate that the near-surface temperature trend pattern is strongly influenced by local sea-ice trends (Jansen et al. 2020;Dai et al. 2019;Boeke and Taylor 2018). For example, the switch of the Barents Sea from a relatively high winter tas trend subregion in the historical period to a relatively low trend subregion under RCP8.5 undoubtedly reflects the early loss of winter sea ice there in many model simulations (Koenigk et al. 2020). This is also reflected in the relatively high tas trend ratios for subregion 13 in Fig. 11d, likely due to the fact that the Barents Sea ice survives farther into the RCP4.5 simulation due to the weaker forcing. It also likely explains why the future RCP8.5 tas trend pattern presented here is markedly different from the pattern obtained by Koenigk et al. (2015) from a comparison of the RCP8.5 tas changes between the periods 2080-2099 and 1980-1999, for the RCA4 subset of these simulations. In that case, the tas changes are very strong over the Barents Sea because of the large increases around the turn of the century that are not included in our "future" trend calculations. This illustrates a limitation of using linear trends to summarize transient climate change and the sensitivity of trend calculations to the choice of endpoints. It should be noted, however, that the CMIP5 historical simulations tend to underestimate the observed Barents Sea ice decline, likely due to the observed Atlantic heat transport across the Barents Sea Opening being on the high end of the regional natural variability (Li et al. 2016). Hence, the simulated timing of the transition from relatively high to relatively low trend may be different from reality.
In cases where the characteristic trend anomaly pattern, noted above, is particularly clear, such as SON tas and JJA sfcWind, the pattern has a stronger amplitude for some model configurations than others. In these cases, and perhaps to a lesser degree cases that exhibit the anomaly pattern less clearly, the subregional trends can be roughly thought of as the sum of the pan-Arctic-Ocean trend, which has a large inter-model spread, and a characteristic anomaly pattern, the amplitude of which varies between model configurations. In fact, for SON tas, as shown in Fig. 8d, most of the raw subregional trends can be reconstructed to within 80% this way, using the weighted mean anomaly pattern with scaling factors ranging from 0.4 to 1.5. While there is some tendency for models with larger anomaly-pattern scaling factors to also have larger Arctic-Ocean means, scaling the raw weighted-mean subregional trends provides a noticeably poorer approximation.
The decomposition of the CMIP5 RCP8.5 surface temperature signal into partial temperature responses, presented by Boeke and Taylor (2018), gives a good framework for discussing this conceptual division of the simulated temperature response into a pan-Arctic component and subregional pattern component, each of which vary in strength between models. While a direct comparison with our results is not possible due to the differing time intervals considered, their focus on surface temperature rather than tas and different seasonal definitions, their results inform a qualitative interpretation of the trend patterns presented here. Consistent with other idealized GCM-based Arctic amplification studies, such as Stuecker et al. (2018) and Pithan and Mauritsen (2014), their annual-average CMIP5 multi-model mean of the partial temperature contributions show a dominant, widespread, positive contribution from clear-sky long-wave radiative effects (LWCS) which include air temperature, moisture and greenhouse gas effects: The spatial pattern of the annual average LWCS is qualitatively similar to the annual average warming pattern. Additionally, all models produce a very similar relationship between the local LWCS contribution and total winter warming. Surface albedo feedbacks (SAF) have the next largest positive contribution, followed by cloud radiative effects. Short-wave clear-sky (SWCS) effects have a small, relatively uniform negative contribution while turbulent heat flux (HFLUX) and ocean heat storage (HSTOR) show strong subregional variation. The contribution with the greatest absolute and smallest relative inter-model standard deviation is LWCS, while HSTOR and HFLUX have a disproportionately large inter-model spread in the sea-ice retreat regions, compared to the sea-ice variation. Boeke and Taylor (2018) report that models that transfer more energy absorbed in summer (primarily through SAF) to fall (primarily through HFLUX) exhibit more Arctic amplification and moreover, those with more effective redistribution of the stored heat warm more overall. This supports the idea that processes determining the temperature trend anomaly pattern (the spatial pattern of LWCS and its relationship to warming and, to a lesser degree, the SAF pattern) are relatively consistent between models, but there is a large inter-model variation in heat dispersal and amplification of this pattern. In fact, Boeke and Taylor (2018) show by regression that, for the CMIP5 models, much of the widespread LWCS heating is related to the HFLUX in the Beaufort/Kara and, to a much lesser extent, Barents/Chukchi Seas. Interestingly, exceptions to this are the Bering Sea and Baffin Bay, which are influenced very little by these major winter ice-retreat regions. The fact that these subregions are relatively unaffected by the dominant cause of inter-model spread in the CMIP5 simulations may explain why they are the subregions in our study that exhibit significantly greater DJF tas trend inter-model spread when the pan-Arctic-Ocean mean is removed. This is not surprising, as both areas are strongly influenced by lower latitudes.
Some qualitative features of the seasonal differences in the tas trend patterns, such as the small JJA trends and larger inter-model spread in DJF than SON, can also be discussed in terms of the individual contributions to Arctic warming. The surface albedo feedback is active only in sunlit months but results in relatively small, local increases in summer warming trends due to the high heat capacity of the ocean. The stored energy delays ice formation in fall and winter, causing strong local tas trends due to increased heat flux from the open water and more widespread positive trends due to temperature and cloud feedbacks, which, in turn, reduce ice formation (e.g., Bintanja and Van der Linden 2013;Boeke and Taylor 2018). This is consistent with our result showing the strongly proportional relationship between summer sea-ice trend anomaly pattern and fall and winter temperature anomaly patterns. Boeke and Taylor (2018) find that the inter-model standard deviation of both the RCP8.5 sea-ice concentration change and surface temperature change in the ice-retreat regions is much larger in January and February than in October, November and December in their CMIP5 analysis. They also show particularly large standard deviations for HFLUX and HSTOR in the ice-retreat regions, suggesting that this is a likely source of the large winter inter-model spread. It would be interesting to see how much of this comes from greater winter natural variability and how much from systematic differences between models.
There is only a hint of this characteristic RCP8.5 trend anomaly pattern in DJF pr, which could be related to the tas pattern via the anticorrelation between clouds and sea ice (for example, Leibowicz et al. 2012). It would be interesting to investigate possible feedbacks in fully-coupled highresolution models, which would benefit from the ability to resolve more realistic physical features than the CMIP5 GCMs and include bi-directional surface-air interactions. Figure 12 shows the DJF and JJA sea-ice concentration (sic) trends and trend anomalies for the driving models. Clearly the JJA trend anomaly pattern is very similar to the reverse of the characteristic anomaly pattern noted above. In fact, as shown in Fig. 13, there is a very strong linear relationship between the JJA driving-model sic trends and the corresponding CORDEX SON tas trends and, to a lesser degree, the DJF tas trends. The relationship between SON tas trend anomalies and JJA sea-ice concentration trend anomalies, illustrated in Fig. 13a, shows remarkable consistency across all the simulations with best-fit proportionality constants close to -0.1 K in all cases, with 68% to 96% of the SON tas trend variance explained by this linear relationship and a rather astounding ≥ 90% for all the ESM-LR-and EC-Earth_r12-driven simulations. The amplitude of the trend anomaly patterns related in Fig. 13a are larger for EC-Earth-driven simulations than the others. Figure 13b shows a wider range of proportionality constants (− 0.06 to − 0.34 K) between the DJF tas trends and JJA sic trends, though they are very similar within driving-model families. Less variance is explained than in the SON case (41-85%), but still more than half for most simulations and at least 80% for those driven by CanESM2.
A simple linear model of the coupled driving GCM temperature and sea-ice trends does a very good job of explaining this relationship. We can write the tas seasonal trend, for  Figure S5a shows the excellent agreement between the actual and best-fit values based on Eq. 1, as well as the fit parameters and r 2 value and Figure S5b shows the residuals, which are small and appear randomly distributed with respect to the trend values. The fact that the best-fit value of SON is negative for NorESM1 might indicate that most of the heat absorbed in summer is stored below the surface melt-water during SON in that model. Similarly, we can represent the SON sic trends as: where D SON represents the spatially uniform component of the stored heat redistribution, SON represents the strength of the first order dependence on remaining heat absorbed by the ocean in summer and the reduced-ice starting point, SON represents the effect of SON atmospheric heating and r SON i is a residual. Figure S6 shows the modelled trends vs fit based on Eq. 2 and corresponding residuals, again, showing a good 1:1 relationship and small, random residuals. Ignoring the residuals, Eqs. 1 and 2 can be combined to give the basis of approximate linear relationship observed in the CORDEX trend anomalies: . The results of a fit to Eq. (3) are shown in Figure S7. As a consistency check, the values of A calculated from the fit parameters of Eqs. (1) and (2) are the same as those calculated by fitting Eq. (3) directly for all models and the values of C are within a factor of ~ 2 except for Nor-ESM1-M. Since the CORDEX model near-surface temperature trends are very close to those of their driving models, this likely explains the observed linearity there as well. These results suggest that many of the complex interactions ignored in this picture average to a relatively small effect on the time and length scales considered here. It might also open up the possibility of characterizing models by the strength of these parameters as a step towards understanding differences in Arctic amplification. A linear relationship with JJA seaice concentration trend anomalies is retained, to varying degrees in the DJF tas trend anomalies (Fig. 13b). For the CanESM2-driven simulations, the proportionality coefficient is -0.34, with 80% to 85% of the variance explained, while for those based on the other GCMs, the coefficient ranges from approximately -0.1 to -0.2 with 40% to 70% of the variance explained by the proportional relationship, suggesting that perhaps the additional summer ocean heating remains more localized for longer in CanESM2 than in the other GCMs, and continues to dominate the near-surface temperature trend for longer. This highlights the importance of local feedbacks in ice-retreat regions. This simple model can be generalized to approximate DJF tas trends as a linear function of JJA sic trends:  1) and (2) to remove the dependence on the SON trends. The scatterplots of DJF tas trends vs the best fits to the JJA sic trends and residuals for the driving models is given in Figure S8, showing impressive linearity. Hence, the trends are strongest during the most rapid phase of sea-ice decline, which varies by region/latitude, which limits the use of linear trends to summarize transient climate change, especially on large spatial scales. A more detailed investigation of these relationships based on cleaner sunlit and dark seasons may give further insight into the relative strengths of the various feedbacks in the GCMs. Given the similarity of the sfcWind trend anomaly patterns to the inverse of that of JJA sea-ice in the COR-DEX models, it is not surprising that there is also a strong apparent inverse relationship between JJA and SON (not shown for CORDEX) sfcWind seasonal trends and JJA sic trends in the CODEX models and driving GCMs. This can be, possibly more tenuously, explained in a similar manner as for tas by assuming that the sfcWind trend is most strongly influenced by the current season's sic trend, since the atmospheric stability and surface roughness are very different over ice than open water. For example, we can write: where W SON is the uniform component of the "atmospheric effect" and j is the residual. Solving Eqs. 1 and 2 for the SON sic trend, assuming the residual is small and substituting into Eq. 5 gives: where Ω = SON SON + SON SON ∕ 1 − SON SON and P � = P SON SON SON ∕ 1 − SON SON + W SON . The results of a fit to Eq. 6 for the three driving models with available sfcWind, are shown in Fig S9. In fact, there is a strong linear relationship between the JJA and SON sic trends (r 2 between 0.72 and 0.97) so it is not surprising that the JJA sic trend pattern is similar to that of SON sfcWind.
A similarly-strong inverse relationship between annual surface wind speed and sea-ice concentration trends for the CMIP5 models, attributed mainly to changes in surface roughness and reduced atmospheric stability, has been reported in Vavrus and Alkama (2021) and their maps of multimodel-mean sfcWind trends are qualitatively similar to our subregional trend pattern, suggesting that our sfcWind pattern is not largely due to the subset of CMIP5 models chosen to drive the CORDEX simulations. Also, Alkama et al. (2020) found, using reanalyses and CMIP5 models, that there is a large-scale feedback of increasing surface winds on sea-ice loss. Likewise, the relationship between SON sfcWind and JJA sic patterns, found here, could could mainly indicate that the near-surface winds are responding to the local sea-ice trends, but it is also possible that they participate in modulating the subregional trend pattern. For example, there could be a positive feedback loop involving stronger summer surface wind trends allowing greater increases in ocean heat storage (Jackson et al. 2010;Timmermans et al. 2018) resulting in stronger negative trends in ice formation and thus positive trends in surface winds (Jakobson et al. 2019;Mioduszewski et al. 2018) and atmospheric heating in fall and winter, followed by further-reduced summer sea ice and possibly strengthened winds (Stegall and Zhang 2012). Of course, in this study, any full feedback loop involving sea-ice can only be represented in the driving GCM, and some air/sea/ice processes are not accurately represented at that resolution. The downscaling models would just be responding to any larger-scale feedbacks in the GCM.

Implications for marine ecosystems
The projected trends in the atmosphere have multi-faceted impacts on the ocean, sea-ice and associated marine ecosystems in the various regions. Since one of the key motivations for this assessment is to evaluate driving models for basin-scale ocean ecosystem models, we highlight some key examples here: 1. Light transmission through the snow-ice system impacts both ice algae phenology and under-ice primary production (e.g., Leu et al. 2011;Lannuzel et al. 2020;Tedesco et al. 2019) and is predominantly driven by the presence of snow and, in spring, the formation of melt ponds. Changes in snow precipitation and particularly snow to rain conversion will have a major impact on light conditions. This is particularly relevant in spring/early summer. Transitions from snow to rain events on sea-ice during that time also enhance the flushing of ice algae and nutrients from the ice into the ocean. As such it is particularly relevant to recognize the impacts of downscaling on the prsn trends. In addition to precipitation, the timing and extent of melt ponds is also impacted by increased and earlier warming (e.g., Light et al. 2008;Abraham et al. 2015) which, in turn, feeds back on both light transmission and flushing. 2. In regions where sea-ice is retreating and the open water area and time is increasing, enhanced warming and enhanced wind and momentum transfer impact phytoplankton growth, nutrient upwelling and carbon uptake (e.g., Steiner et al. 2014;AMAP 2013AMAP , 2018b. The consistent subregional structure for trends in wind and atmospheric temperatures suggest regional differences. For example, we expect these impacts to be less in sub-Arctic seas, such as the Bering Sea, where the ice retreat is projected to be less pronounced in the future and the atmospheric trends are weaker, but particularly high for inflow shelf Seas where the projected trends are relatively strong (Beaufort and Chukchi). Particularly over the Chukchi Sea, wind patterns are suggested to be directly linked to the Bering Strait inflow which impacts the marine ecosystem of the Chukchi and Beaufort Seas (Serreze et al. 2019) in addition to warming and sea ice retreat. These impacts filter through the trophic cascades with a capacity to modify the entire ecosystem of the region (Huntington et al. 2020). 3. A reduced and more mobile ice cover responds more quickly to wind patterns and intensifies the impacts of storms and high winds on wave patterns and coastal flooding. This then affects whale movement (e.g. beluga, bowhead, narwhal), such as migration along the ice edge, areas of congregation and predator avoidance (e.g., Loseto et al. 2018;Scharffenberg et al. 2019;Mathews et al. 2020). In regions where winds show increasing trends in addition to the ice retreat these impacts will be intensified. 4. Wind affects sea-ice drift and deformation as well as upper-ocean mixing and upwelling which directly impacts associated biogeochemical processes. Hence, underestimations in wind speedchanges over the Arctic Ocean under climate change in low-resolution general circulation models can lead to inadequate representation of changes to these processes and the ecosystem as a whole. Often computational resources limit the ability to run ocean ecosystem models with multiple driving models. In such cases, the presented range and uncertainty estimates for atmospheric trends are essential in assessing the uncertainty in corresponding simulations of ecosystem responses.

Caveats
One large source of uncertainty in this analysis is the relative simplicity of the sea ice parametrizations and low resolution of the driving GCMs. The use of prescribed boundary conditions from GCMs, especially for the sea ice and sea surface, is a significant limitation of the CORDEX simulations. Fully-coupled downscaling or high-resolution global coupled models with sophisticated sea ice would be necessary to fully address questions related to seaice feedbacks on a more local scale. These capabilities are on the horizon, with CMIP6 coupled-model transient contributions with a nominal 50 km resolution from one model, CNRM-CM6-1-HR (e.g., Voldoire 2019) and several models contributing 100 km resolution global coupled simulations. Detailed analysis of Arctic climate change is also limited by the lack of spatially-dense, long term observations. Future Arctic observations may be much improved using satellite data, such as that of the planned Atmospheric Imaging Mission for Northern Regions (AIM-North) (Nassar et al. 2019). In addition, developments in statistical interpolation techniques promise added value to the limited observational data of the past (e.g., Newman et al. 2015Newman et al. , 2019.

Conclusions
The Coordinated Regional Downscaling Experiment (COR-DEX) includes a particularly suitable suite of simulations for investigating detailed atmospheric trends in the Arctic. Nine contributions include historical and RCP8.5 scenario atmospheric simulations on the high-resolution ARC-22 and ARC-44 domains. These simulations generally perform well for the recent historical period, 1986-2015 compared to ERA5 and other observations and analyses are mostly only significantly different in subregions 4 and 5, which may be more affected by resolution issues. Observation-based seasonal and annual averaged near-surface air temperature, precipitation and near-surface wind speed trends generally fall within, or near, the model spread, which is comparable to an estimate of the natural variability based on annualaveraged trends in the CanESM2 large ensemble over the Arctic Ocean subregions considered. The analysis provides the following key messages: 1. The RCP8.5 "high-emission" future projections show annual averaged tas trends of between approximately 0.03 and 0.18 K/year from 2016 to 2100, with much stronger trends (and inter-model spread) in winter (0.05 and 0.30 K/year) than summer (0.01 and 0.10 K/year). For subregions representing the Atlantic in and outflow and the Pacific inflow, the subregional trends in all simulations are less than the mean trend of "Arctic-Ocean average" tas, which is a fairly robust indication that these areas experience weaker than average future tas trends. Similarly, the central Arctic and the Canadian Basin, Chukchi and Kara Seas have consistently larger than average tas trends under RCP8.5. The clearest subregional trend pattern can be found in fall. 2. Positive trends in pr and sfcWind are indicated by the majority of the simulations, with mixed results for snowfall trends resulting from the competing effects of warming-derived snow to rain conversion and increased total precipitation. As a result, the decadal-scale behavior of the prsn timeseries are not as linear as for other variables. 3. For tas, pr and sfcWind the RCP4.5 simulations produce approximately 40% weaker trends than RCP8.5, with very similar subregional patterns. The tas and pr multimodel mean ratios are very similar (differing by less than 0.1), indicating that a linear relationship between total Arctic pr and tas change in the RCA4 models 1 3 (Koenigk et al. 2015) also holds for the Arctic Ocean subregions, individually. 4. Comparing simulations that share driving GCMs or downscaling RCMs indicates that the regional tas trends as well as the inter-model spread are generally more strongly influenced by the driving GCM. This grouping suggests that trend uncertainties determined from multi-model and multi-ensemble global modelling efforts might translate to the downscaling model, at least for temperature. This GCM-RCM consistency holds to a limited extend also for the other variables. However, for snowfall and, in some cases, total precipitation, the downscaling RCM has a large effect on local trends and sfcWind trends appear to be systematically higher in the downscaled simulations, regardless of the RCM used. This suggests that downscaling models could have a relatively large effect on trends at a local and seasonal level. This is in addition to downscaling providing more realistic small-scale behavior which may not translate into significant area-averaged trend differences. These are important points to consider when using atmospheric models to provide boundary conditions for ocean and sea-ice models. 5. While the absolute projected trends are often not significantly different between subregions, removing the uncertainty associated with the pan-Arctic trend often exhibits robust relative subregional differences. These qualitative future Arctic subregional trend anomaly patterns show a remarkable consistency among models, scenarios and, in some cases, across seasons and different variables. This suggests that underlying mechanisms, with similar contributions in all models, determine the trend anomaly pattern. 6. There is a strong similarity between the JJA sic subregional anomaly pattern and the SON and DJF tas trend anomaly patterns, as well as the JJA and SON sfcWind anomaly patterns. In fact, the driving-model absolute subregional trend patterns of tas and sfcWind in these seasons can be approximated to a high degree in terms of JJA sic using very simple physically-motivated linear models.
Evaluating the sub-regional variability in trends is important as Arctic climate change and ecosystem impacts cannot be characterized by the Arctic-mean changes alone and it is useful to know which regions are likely to experience consistently higher or lower trends than the still highly-uncertain mean. As indicated above, the qualitative pattern of these trend anomalies is often strikingly similar across variables and between scenarios, with winter tas, both winter and summer sfcWind (and pr, to the extent that the trends are significant) as well as annual snowfall all indicating trend anomalies that are negative or near zero for subregions of the Atlantic sector and in the Bering Sea, and positive in the Siberian and central Arctic Ocean. This pattern is closely tied to the summer sea-ice anomaly trends in the driving GCMs and may be related to feedbacks and connections between variables, such as ice-albedo feedbacks and possibly the anticorrelation of sea-ice with wind speed, and cloud cover, as well as communication between seasons via ocean heat storage and transport. Simple linear models relating subregional summer sea-ice trends and fall and winter tas and sfcWind trends suggests that, to lowest order, much of the complexity of ocean-ice-atmosphere interaction in GCMs has a relatively uniform effect on near-surface atmospheric trends on these time and space scales, and offers a simpler interpretation of differences in Arctic amplification between models.
The trends analyzed here provide key measures for climate change and are associated with many other changes in the Arctic, such as sea-ice retreat, ocean warming and consequent impacts on biogeochemical processes and Arctic marine and coastal ecosystems. These results suggest that Arctic changes already seen over the recent decades will likely continue or accelerate into the future in response to the atmospheric trends.
Code availability Analysis code will be made available upon request. Mask file is in Online Resource 1.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.