Studies have found that ARs are important suppliers of fresh water to many regions around the globe (Dettinger et al., 2011;Rutz & Steenburgh, 2012;Viale et al., 2018). For example, ARs contribute up to 50% of the annual precipitation over California (Dettinger et al., 2011). While ARs can deliver much-needed rainfall in periods of drought, they are also the culprits behind many extreme weather events, such as extreme precipitation (Gao et al., 2016;Lamjiri et al., 2017;Lavers & Villarini, 2013;Ma, Chen, & Guan, 2020;Ma, Norris, & Chen, 2020, extreme wind events (Waliser & Guan, 2017) and flooding (Lavers et al., 2012;Paltan et al., 2017;Ralph et al., 2006;Ruby Leung & Qian, 2009). Besides their influence in the midlatitudes, ARs also play an important role in shaping the polar climate (Baggett et al., 2016;Francis et al., 2020;Nash et al., 2018;Wang et al., 2020;Wille et al., 2019;Yang & Magnusdottir, 2017). For example, Francis et al. (2020) found that ARs influence the formation of Polynyas, which are large and persistent areas of open water and thin ice in thick sea ice-covered areas, of the Southern Ocean. With ARs shifting toward higher latitudes over the Southern Hemisphere in recent decades, it is expected that ARs may exert a stronger influence on the climate of Antarctica (Ma, Chen, & Guan, 2020;Ma, Norris, & Chen, 2020).

Meanwhile, ARs are strongly influenced by the variability of both the large-scale circulation and moisture field (Gao et al., 2015;Payne et al., 2020;Tan et al., 2020;Zavadoff & Kirtman, 2020). Therefore, Arctic sea ice loss may be capable of influencing AR variability in midlatitudes. Yet how ARs are going to respond to Arctic sea ice loss remains unknown. Endorsed by the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016), the Polar Amplification Model Intercomparison Project (PAMIP) provides a guideline for running coordinated multi-model sea ice loss experiments (Smith et al., 2019). PAMIP provides an unprecedented opportunity to study how Arctic sea ice loss influences ARs. Using data from PAMIP, we seek to understand how ARs in the Northern hemisphere winter respond to Arctic sea ice loss, partitioning the response in a dynamical and a thermodynamical contribution. We also investigate uncertainties associated with the AR response to Arctic sea ice loss, using the spread in the multi-model ensemble.

Two sets of experiments from PAMIP are used to investigate how ARs respond to Arctic sea ice loss: "pdSST-piArcSIC" (piArcSIC) and "pdSST-futArcSIC" (futArcSIC). These experiments are atmosphere-only time slice experiments initialized with conditions from April 1st, 2000. Each experiment is run for 14 months with the first 2 months discarded as spin-up. Both experiments are driven by the same present-day sea surface temperature (SST) condition which is based on the 1979-2008 climatology from the Hadley Center Sea Ice and SST data set (Rayner et al., 2003). These two experiments differ only in the Arctic sea ice concentration. piArcSIC is driven by preindustrial sea ice conditions derived from a multi-model ensemble of CMIP5 (Taylor et al., 2012) preindustrial control runs. futArcSIC is forced by sea ice conditions constructed using an ensemble of CMIP5 RCP8.5 scenario, that corresponds to an estimate of Arctic sea ice under  2 C E global warming in these runs. More details on the PAMIP experimental design can be found in Smith et al. (2019). Throughout this manuscript, the response to sea ice loss is calculated by taking the difference between futArcSIC and piArcSIC. The nine models that are included in the PAMIP database at the time of this study are used (See Table S1 in the Supplementary for more Information). Each model experiment consists of 100 ensemble members. Daily specific humidity, zonal wind, and meridional wind at 1,000, 850, 700, and 500 mb from these models and the ECMWF Interim reanalysis (ERA-Interim; Dee et al., 2011) are used to calculate the integrated water vapor transport (IVT) which is then used for AR detection.

We use the Guan and Waliser (2015) AR detection Algorithm. This algorithm is an IVT-based global AR detection algorithm. Many criteria are employed in this algorithm to detect ARs. We will outline some of the key criteria here, and a full description of the algorithm can be found in Guan and Waliser (2015). The algorithm first identifies contiguous regions ("object") of elevated IVT with the regionally and seasonally dependent 85th percentile IVT intensity. The 85th percentile IVT threshold is calculated based on overlapping 5 month windows, so for a given month at a given location, the threshold is derived over all time steps during the five months centered on that month over the entire period covered by the data set. To ensure coherence of the object, the algorithm requires at least more than half of the grid cells within the object to have mean IVT direction deviating less than  45 E within the mean-object IVT direction. Since ARs play a critical role for the poleward moisture transport, the object means poleward IVT is required to be greater than 50 kg/m/s. Lastly, the object should be longer than 2,000 km and with length to width ratio greater than 2. We calculate the IVT threshold separately for each model. The IVT threshold is derived from the piArcSIC experiment and used for both the piArcSIC and futArcSIC experiments for AR detection.

We are interested in the AR response during December, January, and February (DJF) which is the active AR season in the Northern hemisphere. To quantify the dynamical and thermodynamical contributions to the AR response, a scaling method which was previously used in Ma, Chen, and Guan (2020) and Ma, Norris, and Chen (2020) is employed in this study. To obtain the dynamical contribution, the DJF moisture in futArcSIC is scaled so that the DJF mean moisture field in futArcSIC is equal to that of piArcSIC in DJF. More specifically, the DJF specific humidity of futArcSIC at each time step, level, and grid is scaled by the ratio

, where pi E Q and fut E Q are the DJF mean specific humidity for the piArcSIC and futArcSIC, respectively, for the level and grid to which this scaling applies. The scaled moisture field is then combined with the futArcSIC wind field to calculate the scaled IVT which would be subsequently used with the IVT threshold derived from piArcSIC for AR detection. By scaling the futArcSIC DJF moisture field "back" to the piArcSIC DJF moisture field, any change in the AR statistics detected with the scaled IVT field can be treated as the response due to the change in wind field, namely, the dynamical contribution. Similarly, to obtain the thermodynamical contribution, the DJF moisture field in piArcSIC is scaled by the ratio

. The scaled moisture field is then combined with the piArcSIC wind field to calculate the scaled IVT field which will be used for AR detection. Strong moisture/surface air temperature gradients are created at the edge of the strong sea ice loss regions. These strong surface air temperature gradients would modify the wind fields, thus influencing ARs. In conventional dynamic/thermodynamic decomposition, such effects are usually treated as the interaction between the dynamic and thermodynamic. However, the scaling method used in this study is only able to decompose the changes into the dynamic and thermodynamic components, thus not able to account for this kind of covariance between moisture and wind. The linearity of the method can be tested by comparing the sum of the dynamical and thermodynamical effects with the combined effect.

The ensemble mean response of AR frequency is shown in Figure 1a. In this study, AR frequency is defined as the fraction of time a grid point experiences AR condition. Under sea ice loss, the strongest anomalies in AR are found in the climatological AR peak regions, namely, the Pacific and the Atlantic basins. The AR response over both the Pacific and the Atlantic exhibits a tri-pole pattern. Over the Pacific, ARs increase between 30N and 40N and decrease on both sides of this band of enhanced AR frequency. This tri-pole pattern in AR response results in enhanced AR frequency over the Northwest United States (US) and reduced AR frequency over both the Southwest US and Alaska. The drying over Southwest US, including California, is consistent with Cvijanovic et al. (2017) who found that Arctic sea ice loss can induce drying over California in a coupled model. In contrast, over the Atlantic, AR frequency is reduced between 40N and 60N , with enhanced AR frequency on both sides. This pattern leads to enhanced AR frequency over the Iberian Peninsula and drying over Northwest Europe. Note that a significant decrease in AR activity also occurs over central Eurasia.

To better understand the mechanisms behind the multi-model mean AR response, we decompose the total response into a dynamical (Figure 1c) and thermodynamical (Figure 1d) contribution. Over the Pacific, the dynamically induced AR response exhibits a tri-pole pattern similar to the total response, except that the reduced AR frequency in higher latitudes around 60N is more pronounced than the total response shown in Figure 1a. Over the Atlantic, the changes in zonal wind induce a dipolar response in AR frequency, with widespread reduction in AR frequency in high latitudes, and enhanced AR frequency at around 30N . The reduced AR frequency over central Eurasia is dynamically induced. The dynamical response in ARs is consistent with changes in the 850 mb zonal wind (Figure 1e). In response to sea ice loss, zonal wind weakens everywhere around 60N and strengthens at lower latitudes. Over the Pacific, the strengthening of the jet occurs mostly over the jet exit region, resulting in the northeastward extension of the jet, consistent with the findings of Ronalds et al. (2020). Over the Atlantic, the weakening and strengthening of zonal wind occur on the poleward and equatorward sides of the jet, resulting in an equatorward shift of the jet. Such equatorward shift of the jet in response to sea ice loss has also been reported in many sea ice loss studies (e.g., Deser et al., 2010;Peings & Magnusdottir, 2014;Screen et al., 2013;Screen, Bracegirdle, & Simmonds, 2018;Screen, Deser et al., 2018). The thermodynamical contribution is mostly confined to the high latitudes above 50N (Figure 1d). It enhances ARs nearly everywhere, especially over the northern Atlantic and Barents-Kara sea regions. This is consistent with a moistening of the Arctic lower troposphere with sea ice loss (Deser et al., 2015;Figures 1f and S1). Arctic sea ice loss also appears to drive a weak drying of lower latitude areas. However, this drying is only statistically significant off the coast of the Southwest US. Overall, the decomposition suggests that over the regions poleward of around 60N , the dynamical and thermodynamical contributions are opposite in sign and tend to cancel each other, resulting in a relatively small net response there. In the midlatitudes, the response is mostly controlled by the dynamical changes. Our decomposition method nearly reproduces the total response, with only a slight overestimate of the positive response over the Barents-Kara seas, corroborating the validity of the method (Figure 1b).

The atmospheric response to Arctic sea ice loss over the midlatitudes is highly uncertain, possibly due to model deficiency, but mostly due to small signal-to-noise ratio of the induced changes (Cohen et al., 2020;Peings et al., 2021). Despite significant signals emerging in the ensemble mean, these signals are relatively small compared to the noise, as defined from the standard deviation of interannual variability in DJF AR frequency (Figure S2). Over the AR peak regions of the Atlantic and Pacific, the amplitude of the mean response only ranges from 10% to 20% of the noise and reaches a maximum of about 30% over the Barents-Kara Seas. Besides this small signal-to-noise ratio in the response, considerable inter-model spread exists. As shown in Figure S3, while most of the models seem to be able to produce the tri-pole patterns shown in the ensemble mean over both basins, the magnitude and position of the tri-pole patterns vary across models. Further analysis suggests that the inter-model spread is mostly caused by model uncertainty in the zonal wind response, i.e., the dynamical responses (Figures S4 and S6). In contrast, the thermodynamical response is quite robust across models (Figures S5 and S6).

In order to explore reasons for multi-model spread, we apply an empirical orthogonal function (EOF) analysis to obtain the principal spatial patterns of uncertainty in the AR response (Langenbrunner et al., 2015). We first regrid all the model responses on a common grid of 0 9 1 25 . .    , then concatenate the response in each member (900 in total, 9 models, each with 100 members). These EOF analyses are performed across the ensemble member dimension, instead of the conventional time dimension. Since the response uncertainty in each ensemble member is dominated by internal variability, the EOF patterns obtained this way will mostly reflect the uncertainty in the response due to internal variability. We focus on the Pacific and the Atlantic basins. Following Barnes and Polvani (2013), the Pacific and Atlantic are defined as the regions from 0 to 90N and 135E to 125W , from 0 to 90N and 60W to 0 , respectively. To test the robustness of the results, we also varied these boundaries by a few degrees in either direction. The results are not sensitive to small variations in these boundaries.

As shown in Figure 2, the first mode of uncertainty over the Pacific, which explains 25.5% of the uncertainty, describes the inter-member disagreement in the strengthening over the AR core regions and the weakening over the poleward side of the core regions (Figure 2a). The associated uncertainty pattern in zonal wind, which is obtained by regressing the wind anomalies onto the standardized principal component (PC), shows the uncertainty in both the strengthening and meridional shift of the jet (Figure 2c). The second uncertainty mode which explains 13.0% of the uncertainty represents the uncertainty in the narrowing of the AR maximum regions and their northeastward extension. The associated wind pattern mostly reflects the uncertainty in the northeastward extension of the jet. Over the Atlantic, the first uncertainty mode shows a meridional shift of the ARs, which explains about 29.6% of the uncertainty. Accordingly, the associated pattern in zonal wind shows the meridional shift of the jet. The second mode describes the uncertainty in the strengthening and eastward extension of the AR maximum regions. It explains about 14.6% of the uncertainty. The associated wind pattern reflects the uncertainty in the strengthening and eastward extension of the jet (Figure 2h).

We also perform the same EOF analysis on the responses of the 850 mb zonal wind (Figure S7). The results show that the first mode over the Pacific and the first two modes over the Atlantic have spatial patterns that are very similar to the corresponding wind patterns obtained from the EOF analyses on the AR response (Compared Figures S7a, S7e, S7f to Figures 2c, 2g and 2h). Indeed, the correlation between the PCs of the first AR mode and the first wind mode over the Pacific, the first AR mode and the first wind mode over the Atlantic, the second AR mode and the second wind mode over the Atlantic are 0.83, 0.86 and 0.75, respectively. Since the uncertainty in the AR response results from uncertainties in the zonal wind and moisture responses, these high correlations suggest that these three AR uncertainty modes are largely driven by the corresponding zonal wind uncertainty modes. However, the correlation between the second AR mode and the second wind mode over the Pacific is 0.52 which is slightly weaker than the correlation (0.59) between the second AR mode and the third wind mode (Figure S8a). This suggests that both the second and third modes of zonal wind variability may play a role in driving the second uncertainty AR mode over the Pacific.

In the above analyses, we show the uncertainty patterns of the response due to internal variability. However, the response uncertainty can also arise from model differences. To understand the pattern of uncertainty associated with model differences, we perform the same EOF analyses on the ensemble mean AR responses of individual models. Note that Peings et al. (2021) shows that 100 members may be not enough to suppress all the internal variability in these PAMIP experiments. Therefore, the EOF modes obtained this way may still contain a substantial influence of internal variability. Nevertheless, performing EOF analysis on the ensemble mean response of individual models can still give us an estimate for the uncertainty due to model differences. Instead of concatenating the response of individual members together, we concatenate the ensemble mean responses of individual models. Despite having only nine models, the dominant mode of uncertainty over the Pacific and the first two modes of uncertainty over the Atlantic (Figure S9) largely resemble those due to internal variability (Figure 2), though the magnitude is smaller, and the spatial patterns are noisier due to a smaller sample size. Over the Pacific, the first uncertainty mode due to model difference (Figure S9a) mostly describe uncertainty in the strengthening of the AR peak regions. In addition, it also reflects the uncertainty in the northeastward extension. Indeed, its associated wind pattern shows the uncertainty in both the meridional shift and strengthening of the jet, especially over the jet exit regions (Figure S9c). Unlike the second AR mode due to internal variability (Figure 2b) which shows a narrowing and northeastward extension of the AR peak regions, the second mode of model difference (Figure S9b) is quite similar to the first mode, but shifted southeastward, making it project more strongly onto the meridional shift mode. Over the Atlantic, the first and second uncertainty modes due to model differences (Figures S9e and S9f) are the same as those due to internal variability, i.e., they show a meridional shift and a strengthening/eastward extension of ARs, respectively.

Combining the results from the multi-model mean responses (Figure 1) and the uncertainty modes (Fig-

ure 2), one can see that the mean AR response projects onto leading modes of uncertainty over both the Pacific and the Atlantic (Figure 3). For example, the mean AR response over the Pacific projects preferentially onto the second uncertainty mode which represents the northeastward extension of ARs. Over the Atlantic, the mean AR response projects strongly onto the first uncertainty mode which is associated with the meridional shift of ARs. Similarly, the mean 850 mb zonal wind response over the Atlantic also projects onto the uncertainty mode associated with the meridional shift of the winds. The mean zonal wind response over the Pacific, which represents the northeastward extension of the jet, shows no clear preference as to which modes it projects onto. Since these leading modes of uncertainty are mostly driven by internal variability, this implies that they should also represent the leading modes of unforced variability in ARs and the 850 mb zonal winds. Indeed, an EOF analysis of interannual variability in DJF AR frequency and DJF 850 mb zonal winds across all 900 members in the piArcSIC experiment produces leading modes of variability nearly identical to those shown in Figures 2 and S7 (not shown). These results suggest that the response to sea ice loss projects onto natural modes of variability.

Under global warming, generally climate models project a robust poleward shift of the eddy-driven jets over both the Atlantic and Southern Hemisphere (Barnes & Polvani, 2013). Some studies further show that this poleward shift is related to the mean jet state: namely, models with an equatorward-biased climatological jet tend to produce a more poleward shift under global warming (Barnes & Hartmann, 2010;Kidston & Gerber, 2010). Barnes and Hartmann (2010) explain such state dependency by invoking the eddy-mean flow feedback. Since wave breaking can lead to eddy momentum flux divergence and thus weakening of the jet, they argue that the poleward shift of the jet, equivalently, the strengthening of the winds on the poleward side of the jet, is sustained by suppressing wave breaking on the poleward side of the jet. However, due to the sphericity of the earth, the poleward wave breaking is suppressed for a jet that is located nearer to the pole in climatology. Burrows et al. (2017) further examines this relationship in a set of idealized model simulations and confirms the importance of subtropical eddy mixing for the persistence of eddy-driven jet under a range of climate forcings, including tropical and Arctic forcings. Chen et al. (2020) showed that the eddy-driven jet changes under different forcings experience similar circulation feedback in spite of the direction of jet shift. As a result, it is expected that poleward biased jets tend to shift less poleward under climate warming, and equatorward biased jets shift less equatorward under Arctic warming. Given such state dependency in the jet responses to warming, one may ask if the AR response to sea ice loss also exhibits similar state dependency.

Over the Pacific, since the AR response is characterized by a northeastward extension of ARs, we explore whether the degree of extension is related to the climatological AR core longitude. We define AR core longitude as the longitude where the meridionally averaged AR frequency between 20N and 50N is maximum over the Pacific. Indeed, as shown in Figure 4a, models with a climatological AR core located more westward tend to produce a stronger positive response over the AR extension regions defined as the regions between 30N and 50N and  180 E and 230W (blue box in Figure 1e). Similarly, their responses also tend to project more strongly onto the second uncertainty mode over the Pacific which is associated with the northeastward extension (Figure 4c). Over the Atlantic, the degree of AR equatorward shift under sea ice loss seems to depend on the latitude of the climatological AR peak regions, which is defined as the latitude of maximum zonally averaged AR frequency between 60W and 0° (Figure 4b). Models with climatological AR peak regions located over higher latitudes tend to exhibit stronger equatorward shifts and their responses also tend to project more strongly onto the first uncertainty mode which represents the meridional shift of ARs (Figure 4d). Moreover, such relations between the AR response and the AR basic state in the models also exist for the westerly jet (Figure S10), though the relationships are weaker. Over the Pacific, models with a more westward-located climatological jet core tend to simulate a stronger northeastward extension. Over the Atlantic, the westerly jet also tends to shift further equatorward in the models that have a more poleward-located climatological jet.

The results presented in this section suggest a potential emergent constraint for the AR response to sea ice loss: given that the AR response appears to depend on the modeled AR climatology, we may be able to use the observed AR climatology to constrain the spread in modeled AR climatology, and thus narrow the spread in the AR response under sea ice loss (Hall et al., 2019). Following such constraints, observations indicate that the AR peak regions over the Pacific and Atlantic will extend northeastward and shift equatorward under Arctic sea ice loss, respectively (Figure 4). However, the limited number of models used in this study prevent us from reaching a definite conclusion. More models are needed to test the robustness of such relationships between mean state and response. Idealized experiments similar to Burrows et al. (2017) could also be conducted to investigate physical mechanisms behind such emergent constraints.

The response of DJF AR frequency to prescribed +  2 C E global warming-equivalent Arctic sea ice loss has been investigated using a set of atmosphere-only simulations from the PAMIP ensemble. Under Arctic sea ice loss, over the Pacific, it is found that ARs extend northeastward and occur closer to the North American west coast. Over the Atlantic, ARs shift equatorward. The AR response in midlatitudes is mostly governed by the changes in dynamics. Above about 60N , the weakening of the zonal wind tends to reduce AR frequency while the moistening of the atmosphere due to sea ice loss tends to enhance AR frequency. These two effects cancel each other, leading to relatively small negative response over the Pacific sector and positive response over the Atlantic sector.

Inter-model spread in the AR responses is mostly due to the spread in the zonal wind response, with the response in the moisture field being quite consistent across models. We further found that the response uncertainty due to internal variability and model differences resemble modes of internal variability. The uncertainty modes over the Pacific describe uncertainty in the strengthening or weakening of the AR core and northeastward extension of ARs. Over the Atlantic, they represent the uncertainty in the meridional shift and eastward extension of ARs. While the patterns of the leading uncertainty modes due to internal variability resemble those due to model differences, their magnitudes are substantially larger (compare Figures 2 and S9). This suggests that internal variability plays a much larger role than model difference in uncertainties, so that a large ensemble of simulations is required to robustly address the signal-to-noise ratio problem. Analyses of additional 200 members of the SC-WACCM4 experiments (not shown) indicates that the noise associated with internal variability is not completely suppressed in the ensemble mean of the 100 member ensemble, consistent with Peings et al. (2021), and thus more than 100 members are needed to suppress the noise due to internal variability. In addition, the future sea ice loss scenario used here corresponds to 2C warming of the global mean surface air temperature. Different sea ice loss scenarios under stronger warming will also help to strengthen the signal, thus the signal-to-noise ratio. Finally, the noisy patterns shown in the uncertainty modes due to model differences also indicate that more models are needed to better understand the uncertainty due to model differences.

Potential emergent constraints between the AR basic state and the AR response were also explored. Over the Pacific, we found that models with a more westward-located climatological AR core tend to simulate a greater northeastward extension under Arctic sea ice loss. Over the Atlantic, stronger equatorward shifts in response to sea ice loss tend to occur in models that have a more poleward climatological peak of AR activity. Together, these results and the observed AR statistics suggest that Arctic sea ice loss promotes a northeastward extension of ARs over the North Pacific, and an equatorward shift over the North Atlantic. However, the limited number of models used in this study prevent us from making a definite conclusion. Future works using more models with identical sea ice loss forcing are thus needed to test and quantify the robustness of the relationships found in this study. These potential emergent constraints also suggest that it is important to improve the representation of AR climatology in models for AR projection studies. Such improvements can be achieved by increasing model spatial resolution. For example, using the Community Earth System Model version 1.3 (CESM1.3;Chang et al., 2020), shows that CESM1.3 with higher spatial resolution can simulate more realistic AR characteristics compared to the version with lower resolution.

Given the high correlation between the westerly jet and ARs, additional improvements on AR climatology can be obtained by better representation of the jet characteristics, such as mean jet position and strength, in models.

Under climate warming, other factors than Arctic sea ice loss (e.g., tropical SST increases or ocean-atmosphere coupling) can also contribute to the changes in ARs. For example, while the loss of Arctic sea ice shifts ARs equatorward over the Atlantic, the tropical upper tropospheric warming induces an opposite response by shifting the westerly jet poleward (Barnes & Screen, 2015;Butler et al., 2010;Chen et al., 2020;Shaw et al., 2016). Considering the dynamic effect, the robust poleward shift of the westerly jet over the Atlantic under warming suggests that the influences of tropical upper tropospheric warming will dominate over that from the Arctic sea ice loss and shift ARs poleward. However, both Arctic sea ice loss and tropical upper tropospheric warming will also modify the moisture field in their way and thus exert different thermodynamic influences on the AR distribution. Furthermore, the influences of Arctic sea ice loss on ARs and the influences of the tropical upper tropospheric warming on ARs may not be linearly additive, they may interact in some nonlinear ways. Taking these factors together, this "tug-of-war" between Arctic warming and tropical upper tropospheric warming makes AR response to warming uncertain. In this study, we isolate the effects of Arctic sea ice loss on AR changes with two sets of PAMIP experiments which differ only in Arctic sea ice conditions. Future works using PAMIP experiments that include SST changes alone or both Arctic sea ice loss and SST changes with/without ocean-atmosphere coupling are needed to better understand the importance of this "tug-of-war" on ARs.