3.1 Surface energy balance and partial temperature changes

Processes contributing to the surface temperature difference between two experiments are evaluated based on the surface energy balance equation. The basic formulation follows Lu and Cai (2010). The surface energy balance equation for a reference climate is given by

where $S=S^{\mathrm{clr}}+S^{\mathrm{cld}}$ and $F=F^{\mathrm{clr}}+F^{\mathrm{cld}}$ are the downward shortwave (SW) and longwave (LW) radiation at the surface, respectively, with the superscripts, “clr” and “cld” denoting the clear-sky and cloud (total-sky – clear-sky) radiative effects, respectively. The upward LW radiation is given by the Stefan–Boltzmann law, $R=\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}$, where σ is the Stefan–Boltzmann constant and T_s is the surface temperature. The surface emissivity is assumed to be 1. H and L are the net upward sensible and latent heat fluxes, respectively, and Q represents the net downward surface energy flux including the latent heat consumed by snow/ice melting. In the ocean, Q is stored locally or transported. For the difference (Δ) between the two experiments, Eq. (1) becomes

where ΔR_j represents the individual energy terms.

The Stefan–Boltzmann law implies that a larger surface warming (ΔT_s) is required to balance the same amount of energy flux anomaly (ΔR) by emitting LW radiation at a colder background temperature (T_s). Laîné et al. (2016) called this effect the “surface warming sensitivity”, whose importance for the Arctic amplification has been pointed out in multiple studies (Laîné et al., 2009, 2016; Ohmura 1984, 2012; Pithan and Mauritsen, 2014). The warming sensitivity and other energy flux terms may be converted to the same temperature scale (partial surface temperature changes) by

where overbars and dashes represent the global mean and deviations from the global mean (local anomaly), respectively, and

Equation (3) enables the quantification of the effect of a colder winter Arctic requiring more warming to balance the anomalous surface energy flux on the same partial temperature change scale as other components. The left side of Eq. (3) is the simulated surface temperature change. The first, second, and third terms on the right side of Eq. (3) represent local feedbacks evaluated with the global mean warming sensitivity, global mean feedbacks with the local warming sensitivity, and local feedbacks with the local warming sensitivity, respectively. Note that previous studies used the tropical mean in place of the global mean (Laîné et al., 2016; Pithan and Mauritsen, 2014). In Sect. 4.3 and 4.5, each component of the first term is evaluated separately, and the second and third terms are evaluated collectively as the “S-B” effect and “synergy” effect, respectively (Table 3). Accordingly, the surface temperature change formulated by Eqs. (2) and (3) can be written in a more explicit form as

Here, α is computed from the ratio of upward to downward SW radiations at the surface; S, F, H, and L are taken directly from the model output; and Q is computed as a residual of surface heat fluxes (net radiation, sensible heat, and latent heat fluxes). $\sum _j\mathrm{\Delta }{R}_j$ is computed by summing changes in all surface energy flux terms after they are either averaged globally for overbars or the global average is subtracted for dashes. We use the average of T_s from the paired experiments (PI and MH, or HIST and RCP4.5) to calculate $\partial R/\partial T_{\mathrm{s}}$. Although using the average of two experiments or a single experiment for this term has little impact on the results of the current study, we found that the average provided better agreement between the two sides of Eq. (5) for larger perturbations such as a quadrupling of the CO₂ experiment. The diagnosis is made for each grid point and each month. All models are used in this analysis. We note that direct comparisons between different forcing simulations are possible as there is no change in the land–sea mask among the simulations.

Table 3A list of the energy flux terms used in Figs. 5 and 12. Row no. 1 represents the strength of the global mean feedback calculated with local warming sensitivity. Row nos. 2–10 represent the strength of local feedback calculated with global mean warming sensitivity.

Download Print Version | Download XLSX

3.2 Interpretation of surface temperature change at partially ice-covered ocean grid points

The surface temperature archived in the CMIP5 dataset represents the grid-mean “skin” temperature. At the fractionally ice-covered ocean grid points, this variable is a mixture of the sea surface temperature (SST) and ice surface temperature. We assume that the surface temperature T_s at each grid point is reconstructed by

where T_o and T_i are the SST and ice surface temperature, respectively, and A is the ice concentration. The factors contributing to the surface temperature difference for the paired experiments are then diagnosed by

The first and second terms on the right side represent the effect of SST and ice surface temperature changes, respectively. The last term on the right side represents the effect of the ice concentration change, which is weighted by the surface temperature difference between ice and water: the reduction of sea ice cover (ΔA<0) and the exposure of the warmer ocean surface to the atmosphere ($T_{\mathrm{i}}-T_{\mathrm{o}}<\mathrm{0}$) lead to an increase in the grid-mean surface temperature (ΔT_s). In the current analysis, T_o, T_i, and A are obtained from the average of paired experiments. We use T_o in place of T_i for ice-free ocean grids. Only five models (bcc-csm-1, CCSM4, CNRM-CM5, IPSL-CM5A-LR, and MRI-CGCM3) are used for this analysis due to the availability of the required variables, and the consistency of the analysis is verified by agreement between the left and right sides of Eq. (7). The diagnosis is made for each grid point and each month.

3.3 Factors responsible for the model spread

The fractional contribution of individual partial surface temperature changes (or feedbacks in other words) to the inter-model spread of the simulated surface temperature change is given by

where V_j is the fractional contribution and ΔT is the surface temperature change (the subscript “s” in ΔT_s is omitted here). The subscripts j and k denote indices for feedbacks (jth feedback) and models (kth model out of n models), respectively. The overbars denote the average over the models. σ is the inter-model standard deviation of the total surface temperature change. The numerator represents the product of the model spread for each feedback and the model spread for the total feedback, while the denominator represents the ensemble variance of the total feedback. Here, the key points are that (1) V_j accounts for 100 % of the surface temperature change when summed over the feedbacks; (2) positive V_j means that the jth feedback amplifies the model spread, while negative V_j means that it suppresses the model spread. We note that the same formula was used in Yoshimori et al. (2011) and the references therein. The statistical significance of the fractional contribution is tested using the Monte Carlo method by randomly shuffling the model index (k) 10⁵ times. The null hypothesis is that the V_j values neither amplify nor suppress the model spread. When the original V_j is outside the range of the 5−95th percentiles of V_j resulting from the shuffling, it is considered significant. The diagnosis is made separately for ocean and land averages in the Arctic region (north of 60^∘ N). All models are used for this analysis.

Figure 1Multi-model mean (all 10 models listed in Table 2) annual mean surface air temperature response (^∘C): (a) ΔMH and (b) ΔRCP4.5.

Download

4.1 Simulated surface air temperature response

Figure 1 shows the ensemble mean of the annual mean SAT response for ΔMH and ΔRCP4.5. In both cases, the warming in the polar regions is larger than for the rest of the world, particularly in the Arctic. The Arctic-mean response is 0.4 and 3.9 ^∘C for ΔMH and ΔRCP4.5, respectively, whereas the global mean response is −0.2 and 1.9 ^∘C for ΔMH and ΔRCP4.5, respectively (see Table 2 for individual models). This feature reflects the so-called Arctic warming amplification in ΔRCP4.5. The warming at high latitudes and cooling at low latitudes in ΔMH are consistent with the annual mean insolation anomaly caused by the obliquity difference. From this figure, it is unclear whether the Arctic warming in ΔMH is due to forcing and/or feedbacks.

Figure 2Seasonal progress of the zonal mean effective radiative forcing, effective radiative forcing (ERF) (a, b; W m⁻²), and surface air temperature change (c, d, ^∘C): (a, c) ΔMH and (b, d) ΔRCP4.5. The ERF for ΔRCP4.5 is drawn using the data from Yoshimori et al. (2018), and it is computed in the current study for ΔMH. Both ERFs are constructed with a single model, MIROC4m-AGCM (Yoshimori et al., 2018). The surface air temperature changes are the means of all 10 models listed in Table 2.

Download

Figure 2a and b show the seasonal progress of the effective radiative forcing (ERF) for ΔMH and ΔRCP4.5, respectively. The ERF is the top-of-the-atmosphere (TOA) radiation change induced by the forcing constituents and is computed here using the atmospheric GCM (MIROC4m) of Yoshimori et al. (2018), with prescribed climatological SST and sea ice distribution. The ERF for ΔMH was computed by applying the PI and MH insolation to the atmospheric GCM (AGCM) separately, with other boundary conditions held fixed. The TOA net radiation in the MH was averaged for 20 years after a 10-year spin-up and the difference from the PI was taken as ΔMH ERF. The ERF for ΔRCP4.5 was drawn using the data from Yoshimori et al. (2018) in which the time-varying historical and RCP4.5 forcings were applied continuously to the AGCM, with other boundary conditions held fixed. The three-ensemble-member mean of the differences between the 2080–2099 and 1980–1999 averages was taken as ΔRCP4.5 ERF. While this so-called Hansen-style method (Flato et al., 2013; Hansen et al., 2005) is one of the standard procedures for calculating future scenario forcing, e.g., ΔRCP4.5, it is uncommon in paleoclimate applications. With this method, the ERF includes both rapid stratospheric and tropospheric adjustments as well as the land surface response to the instantaneous radiative forcing. Although the land surface response should not be considered as a forcing, we present the ERF to facilitate a consistent comparison between different perturbation experiments. As a supplementary reference, another measure of radiative forcing evaluated by ΔS(1−α_p) is presented for ΔMH in Fig. S1. Here, ΔS is the insolation anomaly and α_p is the pre-industrial planetary albedo. The ΔMH forcing patterns in both Figs. 2a and S1 in the Supplement are qualitatively similar to the familiar insolation anomaly ΔS (e.g., Hewitt and Mitchell, 1996; Ohgaito and Abe-Ouchi, 2007): an increase and a decrease in summer and autumn, respectively, in the Northern Hemisphere, and an increase and a decrease in spring and summer, respectively, in the Southern Hemisphere. For the Arctic average (> 60^∘ N), the peak positive ERF of about 19.9 W m⁻² occurs in July and the peak negative ERF of about −4.8 W m⁻² occurs in September. The ΔRCP4.5 ERF is, in contrast, spatially and seasonally more homogeneous with an annual mean of about 3.0 W m⁻² for the Arctic region. Figure 2c and d show the ensemble mean of the seasonal progress of SAT changes for ΔMH and ΔRCP4.5, respectively. A common and striking feature is that the maximum Arctic warming occurs in autumn (though the magnitude differs substantially) when the ERF is negative or weakly positive. This result suggests that feedbacks play an important role in shaping the seasonality of the Arctic warming for both ΔMH and ΔRCP4.5. This interpretation is in line with Zhang et al. (2010) for ΔMH and Laîné et al. (2016) for ΔRCP4.5.

Figure 3Seasonal progress of the surface air temperature change (^∘C) in the Arctic (north of 60^∘ N): (a) ΔMH land; (b) ΔMH ocean; (c) ΔRCP4.5 land; and (d) ΔRCP4.5 ocean. Thick black lines show the multi-model mean. Note that the range of vertical axis is different for ΔMH (a, b) and ΔRCP4.5 (c, d).

Download

Figure 3 shows SAT changes over the land and ocean for individual models. The seasonality of the SAT change over land is distinct between ΔMH and ΔRCP4.5, but there are some similarities over the ocean: the warming is modest in summer and largest in autumn. Significantly, the model spread over the ocean is also larger in autumn than in summer. The maximum land warming in summer for ΔMH corresponds to the maximum local insolation anomaly, and it thus may appear that the SAT warming over land is not related to the SAT warming over the ocean. However, there are strong cross-model correlations at the 5 % statistical significance level (Student's two-tailed t test) between the Arctic land and ocean for the October–November–December (OND) mean as well as for the annual mean (0.95 for OND and 0.94 for the annual mean). The statistically significant cross-model correlations at the 5 % level also exist for ΔRCP4.5 (0.92 for OND and 0.89 for the annual mean). In addition, the inter-model variance of the Arctic-mean SAT anomaly is larger over the ocean than over land. Although the available surface temperature proxy data for the mid-Holocene Arctic are more abundant on land than over the ocean (Bartlein et al., 2011; Sundqvist et al., 2010), it is useful to focus our analysis on the oceanic region, which has a larger response, and to explore which processes are responsible for the model difference there. We note that there is no statistically significant correlation at the 5 % significance level between ΔMH and ΔRCP4.5 for either the OND or annual means (for both the ocean and land in the Arctic).

Figure 4Surface air temperature anomaly (^∘C) for ΔMH from the simulations (shading) and reconstruction (solid circles): (a) annual mean; (b) July; and (c) January. The reconstruction data are taken from Sundqvist et al. (2010). The mean of all 10 models listed in Table 2 was used.

Download

4.2 Comparison with proxy data

Figure 4 shows the ensemble mean of the simulated ΔMH annual mean, July, and January SAT anomalies superimposed with the reconstructed SAT anomaly at proxy sites taken from Sundqvist et al. (2010). We note that a detailed comparison with earlier PMIP1 and PMIP2 simulations was given by Zhang et al. (2010). There is substantial disagreement between the model and the reconstruction: the warming indicated by the reconstruction is not captured by the model mean in January as well as in the annual mean. The discrepancies are on the order of a few degrees. Although better agreement is seen in July, the simulated warming is overestimated at some North American sites. O'Ishi and Abe-Ouchi (2011) reported that the model–data discrepancy improved substantially when the interaction between the MH climate change and vegetation distribution change is included in one model although the improvement is somewhat limited in other models (Zhang et al., 2010). Unfortunately, none of the models analyzed in the current study include this dynamic vegetation feedback. Comparisons of the model ensemble mean with Bartlein et al. (2011) for the ΔMH annual mean, warmest month, and coldest month are shown in Fig. S2. We note that a more comprehensive comparison with PMIP2 and PMIP3 simulations was presented in Harrison et al. (2014). Again, the model–data discrepancy is large, although the qualitative tendencies of the warming in parts of Scandinavia appear in both. While these limitations need to be kept in mind, they do not reduce the significance of the following results on the understanding of the Arctic warming process. As stated in the introduction, it is not the purpose here to derive a specific emergent constraint using these proxy data, as such a study requires a rigorous statistical approach in parallel to the mechanism understanding and appropriate proxy searches, and is beyond the scope of this article.

Figure 5Simulated and diagnosed surface temperature changes (^∘C) for the ocean (north of 60^∘ N): (a) ΔMH and (b) ΔRCP4.5. The solid black polygonal lines denote simulated changes and dashed blue polygonal lines denote the sum of the diagnosed partial changes; the two lines are superimposed. The graphs represent the means of all 10 models listed in Table 2. See Table 3 for the interpretation of each component.

Download

4.3 Partial temperature changes

Figure 5 shows the contribution of individual energy flux components to the surface temperature change (partial T_s changes) over the ocean in the Arctic diagnosed by the feedback analysis described in Sect. 3.1. As expected, the simulated T_s changes (solid black polygonal lines) are reproduced by the sum of the individual contributions (dashed blue polygonal lines), indicating that the decomposition is useful.

In spring (March–April–May), the total surface temperature change is negative for the case of ΔMH, whereas it is positive for ΔRCP4.5. Therefore, there is no analogy in the response between the two cases. While the synergy effect of local Arctic feedbacks and local warming sensitivity (synergy) slightly contributes to the warming in both cases, the contributions from the downward clear-sky LW radiation components (clr_lw) have opposite signs between the two cases. The albedo feedback (alb) exhibits a relatively large warming effect for ΔRCP4.5, accompanied by cooling due to the surface effect in late spring (net surface heat flux component, or equivalently ocean heat storage and dynamics components). On the other hand, the surface effect is positive for ΔMH, and is accompanied by anomalous turbulent heat fluxes from the ocean to the atmosphere (evap and sens).

In summer (June–July–August), the total surface temperature change is positive but small for both ΔMH and ΔRCP4.5. The albedo feedback is distinctly positive for both cases. An even larger clear-sky SW radiation component (clr_sw) contributes to the additional warming for the case of ΔMH, which is largely driven by the astronomical forcing but it plays little role in ΔRCP4.5. The increased SW radiation reaching the sea surface through the albedo feedback and/or increased seasonal insolation is counteracted by the increased net surface heat flux component, implying that the extra energy is likely stored in the form of ocean heat content. The net result is a small surface warming in summer. It is a common feature of both ΔMH and ΔRCP4.5 that the SW cloud radiative effect (cld_sw) weakens warming by the albedo feedback. This canceling role of clouds in the warm season is consistent with previous studies using future climate projections (Crook et al., 2011; Laîné et al., 2016; Lu and Cai 2009). In both cases, the downward clear-sky LW radiation component plays a substantial role in warming the surface (except for ΔMH in June).

From September to January, the total surface temperature change is larger than in other seasons for both ΔMH and ΔRCP4.5. From September to November, the clear-sky SW radiation component associated with the astronomical forcing contributes to the surface cooling for ΔMH, which is absent for ΔRCP4.5. From October to January for both ΔMH and ΔRCP4.5, the positive surface effect is counteracted by the negative surface turbulent flux components, indicating that the heat is released from the ocean to the atmosphere in the form of latent and sensible heat fluxes. It is, however, unclear how the heat release to the atmosphere leads to the surface warming (or, more precisely, grid-mean skin temperature rise). This point is discussed in the next subsection in detail. It is a common feature of both ΔMH and ΔRCP4.5 that the LW cloud radiative effect (cld_lw) helps warming by the surface effect. This amplifying role of clouds in the cold season is consistent with previous studies using future climate projections (Laîné et al., 2016; Yoshimori et al., 2014). The general increase of cloud cover in autumn to winter for both ΔMH and ΔRCP4.5 is consistent with the enhanced greenhouse effect of clouds (Fig. 6a and c).

Figure 6Seasonal progress of the total cloud fraction change (%) over the ocean (north of 60^∘ N): (a) ΔMH ensemble mean; (b) ΔMH ensemble standard deviation; (c) ΔRCP4.5 ensemble mean; and (d) ΔRCP4.5 ensemble standard deviation. All 10 models listed in Table 2 are used.

Download

Throughout the year, the downward clear-sky LW radiation component exhibits a large contribution and follows the shape of the seasonal progress of the total response for ΔRCP4.5. This component is, however, not large in winter (and June) for ΔMH. This term includes the effect of air temperature and specific humidity changes (and also the radiative forcing of greenhouse gases for the case of ΔRCP4.5), and is qualitatively consistent with changes in both variables (Figs. 7a, c, 8a, c). Obtaining a clear physical interpretation of its role in the surface temperature change is difficult because the primary component of clear-sky LW radiation is emitted from the atmospheric layer near the surface (Ohmura, 2001) where the temperature is tightly coupled with the surface, thus obscuring the causality. Nevertheless, the importance of this component has been reported in previous studies (Pithan and Mauritsen, 2014; Sejas and Cai, 2016). The positive local feedbacks in the cold season with a larger local warming sensitivity make the synergy term an important contributor to the total response for both ΔMH and ΔRCP4.5, as found by Laîné et al. (2016) in future climate projections. For completeness, the same analysis for the land surface temperature is shown in Fig. S3.

Figure 7Same as in Fig. 6 but for the air temperature change (^∘C) (north of 60^∘ N). All 10 models listed in Table 2 are used. (Note that the figure appears blocky compared to Fig. 6 due to the use of a different interpolation scheme in the plotting software which was chosen to avoid a technical issue for pressure coordinate, but it is irrelevant to the data.)

Download

Figure 8Same as in Fig. 6 but for the specific humidity change (g kg⁻¹) (north of 60^∘ N). All 10 models listed in Table 2 are used. (Note that the figure appears blocky compared to Fig. 6 due to the use of a different interpolation scheme in the plotting software which was chosen to avoid a technical issue for pressure coordinate, but it is irrelevant to the data.)

Download

4.4 Interpretation of surface temperature change in partially ice-covered ocean grids

Figure 9 shows the surface temperature change (left side of Eq. 7, ΔT_s) and the individual contributions of surface conditions (the individual terms on the right side of Eq. 7). The surface air temperature change (ΔT_a) is also plotted for reference. The seasonal progress of ΔT_a closely follows that of ΔT_s, suggesting the importance of understanding the grid-mean surface temperature change. The surface and surface air temperature changes have maximum values of 3.3 and 2.9 ^∘C, respectively, in October for ΔMH. They have maximum values of 10.2 and 9.3 ^∘C in November for ΔRCP4.5. The figure indicates that the large increase in grid-mean surface temperature during winter is largely due to the ice surface temperature increase when the SST anomaly decreases seasonally through oceanic heat release after its peak value (Fig. 10a and c). The contribution from ice temperature change has a maximum value of 2.2 ^∘C in October for ΔMH and 6.8 ^∘C in November for ΔRCP4.5. The contribution from SST change has a maximum value of 0.7 ^∘C for ΔMH and 2.0 ^∘C for Δ RCP4.5, both in August. The magnitude of the SST anomaly effect on the grid-mean surface temperature change is small, as the SST change itself is small because SST is fixed at the melting point where sea ice is present and due to the large heat capacity of sea water. The reduction of sea ice cover makes an important contribution to the grid-mean surface temperature increase during autumn. Its peak contribution does not, however, coincide with the timing of the maximum ice concentration anomaly (ΔA; Fig. 11a and c), as the effect is weighted by the surface temperature difference between the sea ice and ocean (T_i−T_o). The interpretation of the results of the feedback analysis in the previous section is that the oceanic heat release in the cold season represented by the positive net surface heat flux term in Fig. 5 contributes to the surface air temperature rise and subsequent ice (and grid-mean) surface temperature rise. This diagnosis is simple but reveals a chain of processes whose temporal links are less clear from the conventional analysis on surface energy balance alone.

Figure 9Contribution of the individual components to the surface temperature change (^∘C) over the ocean (north of 60^∘ N): (a) ΔMH and (b) ΔRCP4.5. The surface temperature change is decomposed into the components of the SST change ((1−A)ΔT_o), sea ice concentration change ((T_i−T_o)ΔA), and sea ice surface temperature change (AΔT_i). Simulated surface temperature (ΔT_s) and surface air temperature changes (ΔT_a) are also plotted for reference. Only five models (bcc-csm-1, CCSM4, CNRM-CM5, IPSL-CM5A-LR, and MRI-CGCM3) are used.

Download

Figure 10Same as in Fig. 6 but for the upper ocean temperature change (^∘C) (north of 60^∘ N). The nine models except for FGOALS-g2 listed in Table 2 are used.

Download

4.5 Factors for the inter-model difference in surface temperature changes

Figure 12 shows the fractional contribution of the partial surface temperature changes to the model spread in the total surface temperature changes. The average is taken for ocean areas in the Arctic, and positive or negative values indicate factors increasing or reducing the model differences, respectively. In the following, individual components whose contributions are either small or inconsistent between the ΔMH and ΔRCP4.5 cases are not discussed, after considering the statistical significance.

Figure 11Same as in Fig. 6 but for the sea ice concentration (%) (north of 60^∘ N). All 10 models listed in Table 2 are used.

Download

In spring (Fig. 12a), large contributions to the model spread are made by the albedo feedback (alb) and the downward clear-sky LW radiation component (clr_lw) for both ΔMH and ΔRCP4.5. Each of these factors contributes to more than 50 % of the model spread. LW cloud feedback (cld_lw) and the synergy effect of local Arctic feedbacks and local surface warming sensitivity (synergy) also contribute to the model spread but to a lesser degree. In contrast, the turbulent heat flux components (evap and sens) as well as the cloud SW radiation component (cld_sw) tend to suppress the model spread.

Figure 12Fractional contribution of individual processes to the model spread in the simulated surface temperature change (%) over the ocean (north of 60^∘ N) for ΔMH and ΔRCP4.5: (a) spring (March–April–May); (b) summer (June–July–August); (c) autumn (September–October–November); and (d) winter (December–January–February) means. The sum of the bar graphs in the same color for each plot adds up to 100 %. The hatching indicates the contribution is statistically significant at the 10 % level. All 10 models listed in Table 2 are used. See Table 3 for the interpretation of each component. Note that the vertical scale for panel (b) is 3-fold larger than those for the other plots.

Download

In summer (Fig. 12b), the albedo feedback (alb) exhibits by far the largest (more than 170 %) contribution to the model spread for both ΔMH and ΔRCP4.5. Note that the vertical scale in Fig. 12b is enlarged 3-fold compared to other plots. As in spring, the downward clear-sky LW radiation component also contributes to more than 50 % of the model spread. The surface effect (net surface heat flux component, or equivalently ocean heat storage and dynamics components) substantially suppresses the model spread for ΔMH, but it is insignificant for ΔRCP4.5.

In autumn and winter (Fig. 12c and d), the downward clear-sky LW radiation component, LW cloud feedback, and surface effect contribute to the model spread, whereas the turbulent heat flux components tend to suppress it for both ΔMH and ΔRCP4.5. As the oceanic heat content is reduced in these seasons through latent and sensible heat fluxes, it is understandable that these two terms have the opposite sign to the surface effect, similar to how the albedo feedback and surface effect have opposite signs in summer. The surface effect contributes to more than 40 % of the model spread in autumn and more than 50 % in winter for both ΔMH and ΔRCP4.5. In contrast to spring and summer, the contribution by the albedo feedback is small in autumn and winter.

The downward clear-sky LW radiation consistently exhibits a large positive contribution (more than 50 %) in all seasons for both ΔMH and ΔRCP4.5. The clear-sky LW radiation is often dominant for ΔMH and ΔRCP4.5 even in spring when the ensemble mean shows surface cooling in ΔMH and warming in ΔRCP4.5. It is also one of major contributors to the model spread even in winter when there is little contribution from the clear-sky LW radiation to the ensemble mean response of ΔMH. The large contribution of this term to the model spread is somewhat expected because this radiative flux largely reflects the surface air temperature, which is thermally coupled with the surface temperature as shown in the previous section. This term, however, also includes the effect of water vapor and lapse-rate changes, whose quantitative contributions are not evaluated separately here. The model variances of the air temperature change are concentrated near the surface in non-summer seasons (Fig. 7b and d), and those of the specific humidity change are large in non-spring seasons (Fig. 8b and d). The relative contribution of air temperature and water vapor to the clear-sky LW radiation may thus vary with the season.

It is also important to point out that the LW cloud feedback contributes positively to the model spread in almost all seasons for both ΔMH and ΔRCP4.5. While the inter-model variability in cloud cover peaks in summer for ΔMH and late autumn for ΔRCP4.5 (Fig. 6b and d), the result suggests that the correct representation of LW cloud feedback is important throughout the year. It is important to recognize that the cloud response is not, however, independent of other feedbacks such as sea ice cover, water vapor, lapse rate, large-scale condensation, and convection (e.g., Abe et al., 2016; Yoshimori et al., 2017). It is also important to notice that the synergy term contributes positively to the model spread. As the surface warming sensitivity depends on the background temperature, this result may suggest that the differences in the reference surface temperature, i.e., model bias, have the potential to reduce the simulated model spread. Taken together, attention needs to be paid to the models' representation of surface albedo, turbulent heat fluxes (and thus the atmospheric stratification including inversion), clouds, and temperature bias to reduce the differences in the models' response.

These results suggest that the processes responsible for the model spread may depend on the season. While the albedo feedback shows only a small contribution to the autumn–winter model spread, this result does not mean that the summer albedo feedback is irrelevant to the model spread in autumn–winter, however. As the reduction of sea ice cover is considered to enhance the oceanic heat uptake through the enhanced albedo feedback, and the reduction of sea ice cover is also considered to enhance the oceanic heat release through the reduced thermal insulating effect, a chain of processes is expected. The model variances of the sea ice concentration change are large from late summer to early autumn with peaks in September–October for both ΔMH and ΔRCP4.5 (Fig. 11b and d), and the model variances of the ocean heat content change are also large in late summer to early autumn, although the peaks occur slightly earlier (Fig. 10b and d). These results are not sufficient to prove the existence of inter-seasonal linkage, but they are consistent with its existence. We calculate cross-model correlations between the summer albedo feedback and October–November–December (OND) feedbacks. The correlations of the summer albedo feedback are 0.72 (ΔMH) and 0.60 (ΔRCP4.5) with the OND surface effect, 0.66 (ΔMH) and 0.69 (ΔRCP4.5) with the OND LW cloud feedback, and 0.85 (ΔMH) and 0.87 (ΔRCP4.5) with the OND surface temperature response (i.e., sum of all feedbacks). These values are statistically significant at the 5 % level according to a Student's two-tailed t test. The significant correlations with the surface effect and with the cloud greenhouse effect are consistent with the chain of processes discussed in Sect. 4.4 and in previous studies (e.g., Abe et al., 2016). Therefore, the model spread in the OND surface temperature response is closely related to the summer sea ice distribution, indicating that feedbacks in preceding seasons should not be overlooked. The recent sensitivity experiment with a single model by Park et al. (2018) demonstrates the dominant influence of sea ice albedo feedback on the MH Arctic winter and annual mean warmings. For completeness, the same analysis for the land surface temperature is shown in Fig. S4.