4.1. Excessive ice-free conditions in CMIP6 models

As reported by Crawford and others (Reference Crawford, Stroeve, Smith and Jahn2021), most CMIP6 models simulate unrealistically long ice-free periods in the HBC (Fig. 2b; see also maps in Supplementary Fig. S1). These biases are often large, with the multi-model mean overestimating the ice-free period by 30 d, and one model overestimating by 87 d. Out of 37 models, 27 (comprising 73% of the models) overestimate the ice-free period, one underestimates the ice-free period, and only nine are unbiased. For every sub-region, a majority of models overestimate the length of the ice-free period. Further, the multi-model mean lies beyond the range of internal variability when compared to the passive microwave record for all regions except for Hudson Strait. Bias is especially strong in southern Hudson Bay, for which only one model simulation falls within the range of internal variability (126 ± 14 d), and 31 of 37 (84%) overestimate the ice-free period by at least 28 d.

Figure 2. Comparison of sea-ice phenology in CMIP6 models and observations (1979–2014) for HBC and its sub-regions (Fig. 1b) using (a, c, e) 80% or (b, d, f) 15% SIC as the concentration threshold (as defined in Fig. 1a). The CMIP6 multi-model mean (red dot) is derived from the first member (or ‘replicate’) of each CMIP6 model historical ensemble (red bars). ERA5 (light gray) and PIOMAS (dark gray) are shown for comparison. The maximum internal variability (2σ _max) from three CMIP6 single-model ensembles is used for the error bars around the mean for the passive microwave record (white dot).

On average, the bias in the ice-free period arises both because sea-ice retreat occurs 19 d too early and sea-ice advance occurs 9 d too late (Figs 2d, f). The greater bias for retreat day is driven mostly by Hudson Strait (retreat bias = −17 d and advance day bias = +4 d) and Foxe Basin (−33 and +12 d, respectively). In other sub-regions, the bias is more equally distributed between retreat and advance, with bias in retreat day being larger by a few days. Eastern Hudson Bay is the only sub-region for which more than one model simulates the retreat day as unrealistically too late, and average retreat day bias is lowest in this region at 8 d too early.

If considering the open-water period (SIC <80%) instead of the ice-free period (SIC <15%), the bias in the CMIP6 multi-model mean is greatly reduced – even eliminated in the northern sub-regions. Most notable is how the opening day (SIC falls below 80%) is timed more accurately than the retreat day (SIC falls below 15%). Therefore, whatever the source of bias, it does not cause the melt season to initiate too soon, but rather to progress too quickly (Supplementary Fig. S2). Similarly, the growth season appears to be too short, starting too late but ending closer to a time that matches observations.

Overall, then, any description of bias source(s) for HBC sea ice must explain four general observations:

(1) The ice-free season is unrealistically long across nearly all sub-regions of the HBC.

(2) The retreat day is consistently too early and the advance day is consistently too late, but the bias in retreat day is greater overall.

(3) The transition seasons (melt season in summer, growth season in fall) occur too quickly.

(4) The magnitude of bias exhibits substantial regional variability.

The following results highlight mechanisms that, if accurately simulated, would reduce the overall model bias in sea-ice phenology. Only mechanisms that proved useful for explaining HBC sea-ice bias are described. Some other mechanisms, which did not prove useful, are mentioned in Section 5.

4.2. Temperature is closely related to sea ice

Average near-surface air temperature over the HBC (hereafter ‘air temperature’) in a model is closely related to its sea-ice phenology (Fig. 3). Many models show a bias in both sea-ice phenology and air temperature, including in the season that precedes a given sea-ice event. Additionally, models that are warmer tend to have earlier retreat and later advance, which leads to a longer ice-free period. This intermodel relationship is roughly linear: if model A estimates the annual air temperature over the HBC as 1°C warmer than model B, we can expect model A's ice-free period to be ~2 weeks longer (beta coefficient in Fig. 3a). This is a strong, significant relationship (r ² = 0.79, p ≪ 0.01). For several model runs, the positive bias in ice-free period can be greatly reduced by correcting for temperature bias. For example, by correcting for the warm air temperature bias of ~+3.0°C in AWI-CM-1-1-MR, its bias in ice-free period can be reduced from ~+58 to ~+15 d. However, several of the models that produce realistic ice-free periods do so by having a cold bias. Correcting the −1.9°C air temperature bias for UK-ESM1-0-LL, for example, would make this model a worse match to observed ice-free periods (bias of −7 d to bias of +20 d). Therefore, air temperature bias is a strong predictor of sea-ice bias overall, but it cannot be the only physical factor driving model bias in the ice-free period.

Figure 3. Scatter plots of HBC sea-ice phenology versus average air temperature (1979–2014). Temperature is averaged annually for the ice-free period (a), during the melt period (b) or the ice-covered period (d) for retreat day, and during the growth period (c) or ice-free period (e) for the advance day. All temperature aggregation is for the HBC region. Black dashed boxes represent the range of internal variability (μ _obs ± 2σ _max). The dotted gray line represents the ordinary least-squares regression of each phenology variable against temperature. The slope of that line is printed at the top of each graph, and an asterisk indicates a significant trend (p < 0.05).

Air temperature bias in a model's 1979–2014 average is also a good predictor of its advance day bias for the same period. The regression line for sea-ice advance predicted by regional temperature in either season passes through the range of internal variability around the observations. This means that if the linear regression of sea-ice phenology on same-season air temperature is used for bias correction, the intermodel spread is reduced and the multi-model mean for average HBC advance day falls within the internal variability range of observations (Figs 4e, f). This improvement comes from bias reductions throughout several sub-regions, especially central Hudson Bay.

Figure 4. Temperature-corrected sea-ice phenology (1979–2014). As in Figure 2, but after applying a bias-correction for regional temperature to sea-ice phenology values in each model. The temperature correction is applied using average annual (a, b), May–July (c, d) or November–December (e–f) temperature.

The average temperature in a model is also a good predictor of its average retreat day; however, the regression line does not intersect the range of internal variability around the observations. In other words, almost all model runs that accurately simulate historical regional temperatures in spring simulate an unrealistically early sea-ice retreat (Fig. 3b). (The one exception is MRI-ESM2-0, which is a good match for both.) Said another way, model runs that accurately simulate the historical average retreat day almost all have a cold bias in May–July. Therefore, if temperature is used to bias-correct the modeled retreat day, modeled retreat day is still too early for most models (70%; Fig. 4d). For some sub-regions, like Hudson Strait or eastern Hudson Bay, applying a bias correction based on temperature actually yields greater bias that is more consistently toward too early retreat. Therefore, although Figure 3b shows that a model's simulated regional temperature has a significant relationship with its simulated sea-ice retreat day, bias in that retreat day cannot be explained solely by bias in average regional temperature.

4.3. Physical processes linking temperature to sea ice

So far, discussion of the air temperature–sea-ice relationship has been strictly statistical, but how does air temperature bias translate to a bias in the timing of sea-ice retreat and (more strongly) advance? A first-order explanation is that higher air temperatures delay sea-ice formation, slow growth and accelerate sea-ice melt. However, there are several complications.

First, since the onset of sea-ice growth depends more directly on the temperature of the ocean rather than the atmosphere, we might expect that a model's sea-surface temperature is more predictive than surface air temperature of its sea-ice advance day. However, May–October sea-surface temperature (r ² = 0.32) is actually less useful than surface air temperature (r ² = 0.51) at explaining intermodel advance day variance (Supplementary Fig. S4).

Second, the timing of sea-ice retreat is related both to the air temperature during the May–July season and the preceding seasons. Models that are warmer during November–April have thinner sea ice in April (Fig. 5a), and that thinner sea ice takes less time to melt come spring (Fig. 6d). A warmer atmosphere may lead to thinner sea ice in two ways: (a) by delaying initial sea-ice advance and (b) by slowing the rate of thermodynamic congelation growth. Models with a warmer November–April do not exhibit significantly slower thermodynamic ice growth (Fig. 5c), but the timing of sea-ice advance is strongly correlated with April thickness (Fig. 6a). April thickness is also strongly correlated with air temperature during November–December (Supplementary Fig. S5), when sea-ice advance typically occurs. Therefore, delayed sea-ice advance is the key reason why warmer models have thinner sea ice in the HBC.

Figure 5. Scatter plots of sea-ice growth versus air temperature (1979–2014) for the entire HBC. Sea-ice growth is defined by (a) the average April thickness, (b) the average rate of change in sea-ice thickness from November to April and (c) thermodynamic thickness change from November to April. Black dashed boxes and gray shading represent the range of internal variability (μ _obs ± 2σ _max). The dotted black line represents the ordinary least-squares regression of the y variables against temperature. The slope of that line is printed at the top of each graph, and an asterisk indicates a significant trend (p < 0.05).

Figure 6. Scatter plots of sea-ice phenology and April sea-ice thickness (1979–2014) for the entire HBC (a, b). In (c), Pearson's correlations between retreat day and the previous year's advance day are on the y-axis and correlations between retreat day and the subsequent advance day are on the x-axis. For models/datasets lying above the 1:1 line (solid gray), retreat day correlates more strongly with the subsequent advance day than the previous advance day. (d) Average sea-ice thickness on the opening day (SIC = 80%) versus the length of the melt period (retreat day–opening day). Black dashed boxes and gray shading represent the range of internal variability (μ _obs ± 2σ _max). The dotted black line represents the ordinary least-squares regression of the y variables against temperature. The slope of that line is printed at the top of each graph, and an asterisk indicates a significant trend (p < 0.05). PMR, passive microwave record.

Third, the strong correlation between air temperature and the ice-free period is due in part to the presence of positive feedbacks between the length of the ice-free period, ocean heat uptake and air temperature (Fig. 6c). For example, when sea-ice retreat is earlier than normal in summer, the ocean can absorb more heat, which both delays the advance of sea ice in fall (Serreze and others, Reference Serreze, Crawford, Stroeve, Barrett and Woodgate2016; Stroeve and others, Reference Stroeve, Crawford and Stammerjohn2016; Crawford and others, Reference Crawford, Stroeve, Smith and Jahn2021) and increases surface air temperature in winter (Screen and Simmonds, Reference Screen and Simmonds2010; Serreze and Barry, Reference Serreze and Barry2011). As shown above, years with delayed freezing in Hudson Bay also tend to yield thinner sea ice, which melts out earlier the next summer (Tivy and others, Reference Tivy, Howell, Alt, Yackel and Carrieres2011). Although it is good that these feedbacks appear in model simulations (Fig. 6c; Crawford and others, Reference Crawford, Stroeve, Smith and Jahn2021), they do complicate proving cause and effect. A warmer atmosphere leads to less sea ice, but a reduction in sea ice also facilitates a warmer atmosphere.

One line of evidence that helps support the idea that excessive warmth is forcing sea-ice bias in many models is the prevalence of southerly near-surface wind biases – especially during the melt season and ice-free season (93% of models). Hudson Bay is dominated by northwesterly winds that blow off mainland Nunavut (Fig. 7), so southerly wind biases indicate that the models simulate unrealistically weak northwesterlies and/or unrealistically frequent wind reversals. In either case, such biases lead to less cold air advection than portrayed in ERA5. Strong southerly wind bias is especially notable in several models exhibiting strong negative biases in sea-level pressure to the northwest of Hudson Bay, like AWI-ESM-1-1-LR and NESM3 (Fig. 8). However, if negative sea-level pressure anomalies to the north are paired with strong positive sea-level pressure anomalies to the south (e.g. BCC, CESM2 and NorESM2 models), the bias in wind direction is more westerly, and models exhibit no temperature bias or a cold bias. By contrast, models with positive sea-level pressure anomalies to the north of Hudson Bay and negative sea-level pressure anomalies to the south (e.g. IPSL and MIROC models), tend to have more easterly wind bias and high-temperature bias. Finally, the two models dominated by positive sea-level pressure anomalies over Hudson Bay in summer (ACCESS-CM2 and KIOST-ESM) are two of the coldest.

Figure 7. Model bias in 2 m air temperature and 10 m wind vectors during August–October (1979–2014), which roughly corresponds to the ice-free season in the HBC.

Figure 8. Model bias in mean sea-level pressure during August–October (1979–2014), which roughly corresponds to the ice-free season in the HBC. Letters in the upper-left corners indicate whether the model exhibits a cold bias (C), warm bias (W) or no significant surface temperature bias (N) compared to ERA5 during August–October (1979–2014).

Additionally, warm biases of the land surface, especially to the northwest, can also lead to underestimation of cold air advection. Notably, the EC-Earth3 family of models, which have minimal bias in air temperature and surface winds, also simulate accurate advance day (Fig. 3e). Together, these findings suggest that the biases in advance day can be traced back to biases in atmospheric circulation during the ice-free season that leads to less advection of cold air over the HBC. In other words, summer temperature bias is a cause of unrealistically late sea-ice advance, not just a correlated variable.

4.4. Sea ice dynamics and regional differences in sea-ice phenology

Sea-ice phenology in the HBC exhibits contrasting spatial patterns of sea-ice advance and retreat. In fall, sea ice starts forming in Foxe Basin and northwestern Hudson Bay before advancing into the south and finally the east (Fig. 1h). If thermodynamics were the only control on sea ice, we would expect the opposite spatial pattern for sea-ice retreat. Counterintuitively, retreat often begins in the northwest in late June, around the same time as eastern Hudson Bay, and retreat in southern Hudson Bay typically occurs ~3 weeks later (Fig. 1g). Sea-ice retreats later in the south than the northwest because of the cyclonic (counterclockwise) circulation within Hudson Bay during winter that leads to substantial divergence of sea ice out of the northwest and convergence and thickening of sea ice in the south (Fig. 9; Saucier and others, Reference Saucier2004; Kirillov and others, Reference Kirillov2020). Note that although PIOMAS is a reanalysis, its sea-ice retreat days and drift are consistent with observational data except in James Bay (Fig. 2; Etkin, Reference Etkin1991; Kirillov and others, Reference Kirillov2020). The eastern side of Hudson Bay also experiences net convergence of sea ice, but its sea ice typically breaks up earlier than the south, which has been linked to substantial input of relatively warm river discharge in spring (Jones and Anderson, Reference Jones and Anderson1994; Ridenour and others, Reference Ridenour2019). Therefore, another way to assess the accuracy of simulated sea-ice phenology is to examine the differences in sea-ice advance and retreat between three key regions: the northwest, the east and the south (Fig. 10).

Figure 9. Decomposition of ice growth in Hudson Bay during the ice-covered season (November–April; 1979–2014) into (a) total, (b) advective growth, (c) convergent growth and (d) thermodynamic growth. Ice ‘growth’ is measured as the average daily change in effective thickness. Vectors indicate average sea-ice drift. Note that this decomposition requires sea-ice drift vectors, which have NaN values along the coast. Data from PIOMAS.

Figure 10. Regional differences in the timing of (a, b) retreat and (c, d) advance day in the HBC (1979–2014). Positive values indicate that the first sub-region in a given pair experiences later retreat or advance. (e) Sub-region definitions. Symbology matches Figures 2 and 4.

For sea-ice advance, all models correctly simulate the freeze-up starting in the northwest (positive values in Fig. 10c), but most (62%) simulate the lag between advance in the northwest and advance in the south as being too long (average bias = 5 d). The multi-model mean lag in advance between the northwest and east is accurate, and about the same number of models overestimate (24%) and underestimate (22%) the lag time. The timing of sea-ice advance in the east versus the south is less accurate; most models (59%) exhibit significant negative bias, and 32% incorrectly simulate advance occurring in the east before the south (negative values). Applying a bias correction based on regional temperature (as for Fig. 4) eliminates the multi-model mean bias in advance day regional differences. Excessive summer heat is more common in the south than either the northwest or east (Fig. 7), which aligns with the south being the region with the worst bias in advance day (Fig. 2f). Therefore, applying the temperature-based bias correction makes advance day earlier in all three regions, but especially the south.

Models exhibit a wide spread in regional differences of sea-ice retreat, but most CMIP6 models (97%) reproduce the observed pattern of earlier retreat in the northwest and later retreat in the south (positive values in Fig. 10a). By contrast, about half of the models (46%) accurately simulate sea-ice retreat occurring at about the same time in the northwest and east, and fewer (27%) accurately simulate that retreat day occurs ~17 d later in the south than in the east. In fact, 40% of models actually simulate retreat occurring in the south first. In other words, biases in retreat day patterns are worse than biases in advance day patterns. Applying the temperature-based bias corrections does improve estimates of the east versus northwest difference (average bias of +10.8 to −6.4) and east versus south difference (mean bias of +13.5 to +8.1), and only 14% models still exhibit a pattern of retreat occurring in the south before the east. However, accounting for average model temperature is less effective at reducing bias in sea-ice retreat patterns than sea-ice advance patterns. Additionally, we find that the bias in the retreat day difference between the south and northwest worsens if a temperature-based bias correction is applied (−2.4 to −14.2). Models are too warm on average in both the northwest and south (Supplementary Fig. S6), so an increase in model bias after applying the bias correction suggests that the high-temperature bias was counteracting bias in some other mechanism that influences sea-ice retreat, such as sea-ice transport between the northwest and south (i.e. sea-ice dynamics). Biases in sea-ice dynamics may also help explain the bias in retreat timing between the east and south or east and northwest, which is not completely resolved by accounting for temperature.

Looking first at maps of sea-ice circulation, 44% of models simulate a cyclonic circulation that is too weak and 33% simulate a cyclonic circulation that is too strong (Fig. 11). This variability in sea-ice drift leads to differences in sea-ice convergence. Recall that PIOMAS sea-ice growth in the northwest of Hudson Bay is primarily thermodynamic, whereas sea-ice growth in the south and east is primarily from convergence. These spatial patterns are summarized by calculating the difference in convergent sea-ice growth between the three sub-regions and comparing to the regional differences in sea-ice retreat (Fig. 12).

Figure 11. Model difference from PIOMAS in April sea-ice thickness and January–July drift (1979–2014). January–July roughly corresponds to the ice-covered season and melt season in the HBC.

Figure 12. Relationship between regional differences of sea-ice retreat and convergent sea-ice growth (January–July, 1979–2014) for (a–c) the northwest, south and east regions of Hudson Bay (defined in (d)). Black dashed boxes represent the range of internal variability (μ _obs ± 2σ _max). The dotted gray line represents the ordinary least-squares regression of retreat day difference against wind velocity. An asterisk indicates that a coefficient is significant at p < 0.05.

In comparison to PIOMAS, all models similarly show more convergent sea-ice growth in the south than the northwest, but over half (56%) simulate a significantly higher difference (Fig. 12a). All else being equal, an overestimate of convergent growth in the south ought to lead to greater sea-ice volume in the south than the northwest, and therefore overestimation of the lag between retreat day in the northwest and retreat day in the south. However, recall that lag in retreat days between northwest and south in model simulations is actually quite accurate (Fig. 10a). How can this be? The high-temperature bias is also greater in the south than the northwest for most models exhibiting excessive convergent growth in the south (Fig. 7), so it is possible that regional differences in the temperature bias and convergence bias have compensatory effects. This would also explain why bias-correcting sea-ice retreat day based on temperature (but not sea-ice dynamics) leads to greater bias in simulations of the northwest to south lag in sea-ice retreat. However, the lack of a significant relationship for intermodel variance in Figure 12a weakens confidence in this explanation.

By contrast, the difference in sea-ice convergence in the east versus other regions shows both a clear bias and a clear relationship with spatial differences in retreat day (Figs 12b, c). Every model underestimates the difference in convergent growth between the eastern region and either the south or northwest. In fact, about half (52%) show more convergent growth in the south, whereas PIOMAS depicts more convergent growth in the east. Moreover, the more convergent growth that occurs in the east compared to other regions, the later the sea-ice retreats in the east compared to other regions. This is exactly as expected since more convergent growth means greater sea-ice volume.

We should use some caution comparing to PIOMAS because, although PIOMAS assimilates SIC, is does not assimilate sea-ice motion or thickness observations. Therefore, the true behavior of Hudson Bay sea-ice dynamics may lie somewhere between PIOMAS and CMIP6 values. That said, the sea-ice phenology values derived from PIOMAS are a better match to passive microwave observations than most CMIP6 models, and causes of the large differences in the sea-ice convergence in PIOMAS versus CMIP6 are still worth exploring.

4.5. Wind biases help explain biases in sea-ice convergence

The ability of a model to reproduce observed regional differences in convergent sea-ice growth relates, at least in part, to how closely it reproduces observed wind fields (Fig. 13; Supplementary Fig. S6). For example, although most models agree with observations that sea-ice retreats in the northwestern region before the southern region, this spatial difference is larger in models with more strongly westerly winds over the western and northwestern regions (Fig. 13a). Since the average wind direction is northwesterly, more westerly means winds are stronger and more effective at pushing sea ice out of the northwestern region and into the southern region. However, note that there is actually an easterly wind bias for several models for which the convergence in the south is too strong. Convergent growth depends on fields of both sea-ice drift and sea-ice thickness, so wind bias provides an important but incomplete explanation.

Figure 13. Relationship between regional differences in convergent sea-ice growth and regional winds. Comparisons are separated by (a–c) zonal winds and (d–f) meridional winds during the ice-covered and retreat season (January–July) 1979–2014. (g) Region definitions. Black dashed boxes represent the range of internal variability (μ _obs ± 2σ _max). The dotted gray line represents the ordinary least-squares regression of convergence difference against wind velocity. An asterisk indicates that a coefficient is significant at p < 0.05. Major outliers are excluded from several regression line calculations: MIROC-ES2L (d, e) and the CMCC models (a–c).

Considering the difference between the east and northwest regions, we find that both zonal and meridional winds are important (Figs 13b, e). Especially if the outlier CMCC models are excluded from calculations, models with more westerly winds on the west side of Hudson Bay also have a greater lag time between retreat in the northwest and retreat in the east. Models also show a greater difference between retreat in the northwest and the east if the meridional winds on the east side of the bay are less northerly. In other words, the more cyclonic the winds are in a model, the longer lag between retreat in the northwest and retreat in the east. All models but one are consistent with ERA5 in simulating northwesterly winds with cyclonic vorticity across Hudson Bay, and there is no clear bias in zonal winds on the west side of Hudson Bay (Figs 13a, b). However, 91% of models have a southerly wind bias on the west side of Hudson Bay (Fig. 13d). On the east side of Hudson Bay, 68% of models have an easterly wind bias and 55% have a southerly wind bias (Figs 13d, e). The direction of these biases all indicate winds are too weak compared to ERA5. One caution, though, is that the NCEP-NCAR reanalysis disagrees with ERA5 regarding the average meridional wind velocity over Hudson Bay. ERA5 is generally the best reanalysis product for Arctic surface wind (Graham and others, Reference Graham2019a, Reference Graham, Hudson and Maturilli2019b), but the NCEP-NCAR reanalysis is used for PIOMAS. This contributes more uncertainty to our statements about meridional winds.

Finally, although satellites and PIOMAS say the southern region should experience retreat after the eastern region, models show no consistent lag (Fig. 10a). Similarly, models are split on whether the south or east experiences more convergent growth (Fig. 12f) despite PIOMAS showing a clear tendency for more convergent growth in the east. However, the importance of wind bias to these regional differences is unclear because there is no significant relationship for the intermodel variance of surface winds and regional convergent growth (Figs 13c, f). Additionally, an easterly zonal wind bias relative to ERA5 (which 91% of models exhibit) ought to decrease convergent growth in both the south and east by reducing transport from the west and northwest.

In summary, the difference in sea-ice convergence between the northwest and south/east regions in a model is significantly related to zonal wind velocity over the west side of Hudson Bay. This helps explain some of the biases in regional sea-ice retreat, especially between the northwest and east. However, biases in simulation of regional differences between the south and the east cannot be explained by winds.