The ice factory of Hudson Bay: Spatiotemporal variability of the Kivalliq Polynya

Located in the Arctic and subarctic regions of Canada, Hudson Bay is a large, shallow, inland sea that, together with Foxe Basin in the north and James Bay in the south, forms the largest inland sea in the world, covering approximately 804,000 km² (Gagnon and Gough, 2005; Hochheim and Barber, 2010). As an inland sea, Hudson Bay is isolated from the open ocean circulation of the Arctic and Atlantic Oceans (Ingram and Prinsenberg, 1998); thus, variations in the seasonal ice cover are largely influenced by atmospheric forcing rather than the advection of ice or water from other basins (Saucier and Dionne, 1998). Within Hudson Bay, freeze-up typically occurs from late October to November in a northwest to southeast pattern, while retreat begins in May and progresses from the northwest to southeast (Gagnon and Gough, 2005; Hochheim and Barber, 2010; Andrews et al., 2017) with the last remnant ice melting out by early August in either southern or eastern Hudson Bay depending on the prevailing wind direction during winter (Kirillov et al., 2020). The ice cover is comprised of a relatively narrow coastal band of landfast ice, though a vast majority of the ice cover is mobile and circulates cyclonically within the Bay (Landy et al., 2017). This motion gives rise to 2 features: (1) a recurrent polynya in northwestern Hudson Bay that arises due to northwesterly offshore winds pushing the mobile ice away from the landfast ice edge and (2) a west–east gradient in ice thickness due to the formation of new, thin ice within the polynya and the deformation of ice into thicker ridges in eastern Hudson Bay (Landy et al., 2017). The latter feature was first identified by Landy et al. (2017) and was recently explored further by Kirillov et al. (2020); however, the polynya in northwestern Hudson Bay has yet to be studied or characterized in terms of its spatiotemporal variability and contribution to both the ice mass balance and physical oceanography within Hudson Bay. In previous research, this polynya has simply been referred to based on its location in northwestern Hudson Bay (e.g., Saucier et al., 2004; Gunn, 2014; Kirillov et al., 2020); however, as it adjoins the Nunavut administrative region of Kivalliq, we elect to call it the Kivalliq Polynya.

Generally, a polynya is defined as a persistent and recurrent area of open water and/or thin ice that occurs at times and under climatological conditions when one would expect a thicker or more consolidated ice cover (Smith et al., 1990; Tamura and Ohshima, 2011). Polynyas may be further categorized based on their mode of formation: “latent heat” polynyas and “sensible heat” polynyas. A latent heat polynya occurs when the ice within a polynya is continually removed by winds or currents and the heat required to balance that loss to the atmosphere is provided by the latent heat of fusion of the ice that forms in the open water (Smith et al., 1990; Morales Maqueda et al., 2004). In contrast, a sensible heat polynya is maintained through the presence of warm waters at depth that often inhibit the formation of an ice cover through vertical mixing or upwelling (Barber and Massom, 2007), though in some instances the heat flux from depth may not be sufficient to completely prevent ice formation but may slow ice growth and create what has come to be known as an “invisible” polynya (Melling et al., 2015). Although both of these mechanisms often play a role in the formation of a polynya, the Kivalliq Polynya is considered a latent heat polynya (Landy et al., 2017).

Morales Maqueda et al. (2004) outline 4 ways that polynyas play an important role in climate. First, the significant ocean-to-atmosphere heat and moisture loss from within a polynya warms the air column above and modifies mesoscale atmospheric processes. Second, the baroclinic circulation of the ocean near the polynya is affected by the increased salinity of the upper ocean resulting from brine rejection during the formation of frazil ice within the polynya. Third, as a result of vertical oceanic mixing and convection, polynyas affect biogeochemical air–sea fluxes. Fourth, polynyas act as a biological hotspot for marine life, spanning all trophic levels from phytoplankton to marine mammals. Studies have found that biological productivity is greater in a polynya than in adjacent pack ice, resulting in the increased abundance of marine mammals and birds (Stirling, 1997). Given their importance in the climate system, Morales Maqueda et al. (2004) suggest that variability in the size and occurrence of polynyas may be useful as indicators of potential changes in climate.

The Kivalliq Polynya has been shown to play an important role in the climate of Hudson Bay. Modeling studies have demonstrated that the polynya is an area of enhanced ice production (Saucier et al., 2004), with the associated brine rejection increasing salinity in western Hudson Bay (Burt et al., 2016) and driving deep-water formation within the Bay (Stewart and Barber, 2010; Granskog et al., 2011). Polynyas are known to be biologically important (Stirling, 1980): specifically, the Kivalliq Polynya affects the seal population and in turn the Western Hudson Bay polar bear population (Stirling and Derocher, 1993; Stirling, 1997). Further, the polynya presents potentially hazardous thin ice conditions along the landfast ice edge for local Inuit who reside along the northwestern coast of Hudson Bay and who rely on this landfast ice for a hunting platform and transportation network (Babb et al., 2019).

Due to the remoteness of Arctic polynyas, in situ observations are difficult; thus, remote sensing has become the primary method used to monitor their development. Visible and thermal infrared satellite imagery is useful for deriving polynya characteristics at a high resolution (approximate meter scale) but suffers from near ubiquitous cloud cover (Markus and Burns, 1995; Cheng et al., 2017). Thus, to create a continuous, long-term record of polynyas, researchers have made use of passive microwave satellite imagery. This imagery is available at daily global coverage and is not affected by cloud cover, though it has a much coarser resolution (approximate kilometer scale; Markus and Burns, 1995). However, this resolution is adequate for observation of the Kivalliq Polynya which has a spatial scale of thousands of square kilometers.

Historically, sea ice concentration (SIC) has been used to define polynyas within an ice cover (e.g., Barber and Hanesiak, 2004; Cheng et al., 2017; Dale et al., 2017); however, there are 2 disadvantages of using SIC to identify and determine the area of a polynya. First, any moisture at the surface, including melt ponds, may be detected as open water and lead to overestimation of the polynya area. Second, if air temperatures fall below freezing, then new ice may begin to form almost immediately within latent heat polynyas and subsequently lead to the measurement of 100% SIC and a complete thin ice cover (Comiso et al., 1997). To overcome these issues, researchers have instead used the polynya signature simulation method (PSSM; Markus and Burns, 1995) or a thin ice algorithm (e.g., Martin et al., 2004a; Tamura et al., 2007; Paul et al., 2015).

We have opted to use the thin ice algorithm adapted by Martin et al. (2005) from that derived by Martin et al. (2004a); this adaptation uses the 12.5-km resolution Advanced Microwave Scanning Radiometer (AMSR) 36 GHz channel. Although PSSM offers a higher spatial resolution (6.25 km) due to its use of the Special Sensor Microwave/Image (SSM/I) 85 GHz channel, it only classifies polynya ice thickness into discrete categories of open water, thin, and thick ice; in contrast, the thin ice algorithm using the 37 GHz channel identifies areas of thin ice and estimates the thickness, but at a lower spatial resolution (12.5 km). Further, a number of thin ice algorithms have also used AMSR-E brightness temperatures (e.g., Iwamoto et al., 2014; Nihashi et al., 2017; Nakata et al., 2019); indeed, Iwamoto et al. (2014) mapped ice production in polynyas across the Arctic Ocean and found a distribution in northwestern Hudson Bay very similar to that reported in this work. However, these algorithms use the AMSR-E 89 GHz channel which suffers from sensitivity to water vapor, cloud cover, and ice fog and may also misclassify fast ice as thin ice (Nihashi and Ohshima, 2015).

The thin ice algorithm derived by Martin et al. (2004a) is based on the ratio of the SSM/I 25-km resolution 37 GHz vertical and horizontal brightness temperatures and is able to determine a continuous pattern of ice thicknesses, with the 10-cm contour used as the boundary of the polynya. Incorporating wind speed, air temperature, and radiation variable data from the European Centre for Medium-Range Weather Forecasts Reanalysis Fifth Generation (ECMWF ERA5) dataset (C3S, 2017) into an energy balance and conductive heat flux equation, this technique allows for the quantification of heat flux through the ice (under assumption of negligible flux from underlying ocean; Drucker et al., 2003; Martin et al., 2004a) and therefore the daily volume of ice production. Martin et al. (2005) subsequently adapted this thin ice algorithm for use with the 12.5-km resolution AMSR-E 36 GHz channel, applying it to the Chukchi Sea polynya. We have taken the additional step of applying this algorithm to the 12.5-km AMSR-2 36 GHz channel in order to extend the record of the polynya to 16 years. Through individual comparison of AMSR-E and AMSR-2 with SSM/I, the findings of Nihashi et al. (2017) suggest that a slight bias exists between the two satellites; however, as their specifications are almost identical, we have opted to apply the same algorithm to AMSR-2 as AMSR-E (Iwamoto et al., 2014). Within this work, we applied the thin ice algorithm to 16 years of AMSR data (2002–2019) to characterize the spatial and temporal variability of the Kivalliq Polynya and assess its impact on the ice mass balance and marine environment of Hudson Bay.

The dataset used in this research is 16 years in length (2002–2019), comprised of 12.5-km resolution AMSR-E and AMSR-2 gridded brightness temperature data at 36 GHz, provided by the National Snow and Ice Data Center (Cavalieri et al., 2014b; Meier et al., 2018). Following previous polynya studies (e.g., Martin et al., 2004a; Kwok et al., 2007), our analysis focused on ice production and polynya area during the winter season (December–April). From our calculations, ice production in the region in December may be attributed to initial freeze-up as well as to the opening of the polynya later in the month. In April, ice production may be attributed to polynya events, but this month also includes the initial stages of ice breakup in the region. Finally, the analysis ran from 2002 to 2019 with a gap during 2012 between the failure of AMSR-E in October 2011 and the launch of AMSR-2 in May 2012.

The study area was defined based on the work of Gunn (2014) who divided the Hudson Bay region into nine subregions. One of these subregions was the Northwest Shore subregion (Figure 1), which extends from Cape Churchill in the south to the southwestern corner of Southampton Island and westward across the mouth of Roes Welcome Sound. While coastal flaw leads extend beyond these boundaries, areas of open water and thin ice associated with the Kivalliq Polynya are generally contained within these bounds. The gridded AMSR datasets were subset to grid cells contained within this region for further analysis of the polynya.

In this study, the polynya open water area was identified using the thin ice algorithm derived by Martin et al. (2004a). This algorithm uses the ratio of the SSM/I 37 GHz channel vertical and horizontal brightness temperatures, R₃₇,

Following Martin et al. (2005), we adapted this algorithm to use the AMSR 36 GHz vertical and horizontal brightness temperatures at 12.5-km resolution. However, because of the daily and seasonally dependent drift in calibration of the AMSR data relative to the SSM/I data (Martin et al., 2004b, 2005), the AMSR and SSM/I data sets had to be cross-calibrated. We obtained the AMSR 36 GHz and SSM/I 37 GHz brightness temperature data both at 25-km spatial resolution (Cavalieri et al., 2014a; Markus et al., 2018; Meier et al., 2019) and then compared spatially coincident values of R₃₆ and R₃₇ in a daily scatterplot (Figure 2). Based on the daily linear regression of these two datasets, R₃₆ was corrected to R₃₇ with the following equation,

With R₃₆ corrected to R₃₇, the polynya ice thickness (h_T) can be calculated using the equation derived by Martin et al. (2004a),

where a comparison with AVHRR gives α = 230.5, β = −243.6, and γ = 1.008 (Martin et al., 2004a). This curve is valid for R₃₇ < 1.4, or for a minimum thickness of approximately 5 mm, and works best for ice thicknesses in the range of 5 mm to 10 cm (Martin et al., 2005). In their study of the Chukchi Sea Alaskan coast polynya, Martin et al. (2004a) found that approximately 90% of heat loss from the polynya came from ice thinner than 10 cm and thus defined the polynya as an area with ice thinner than 10 cm. Other polynya studies have used a 20-cm threshold (e.g., Paul et al., 2015; Fraser et al., 2019); however, this threshold was based on an ice thickness algorithm derived by Adams et al. (2013) that used MODIS ice surface temperatures in combination with NCEP/DOE Reanalysis II and was considered appropriate for the retrieval of thicknesses up to 20 cm. As we are using the thin ice algorithm derived by Martin et al. (2005), we have opted to use their threshold of 10 cm or less. Using this threshold, the total daily polynya area was then calculated through 16 winters.

Because new ice has little to no snow on the surface and assuming that the ice–air interface is not melting, the heat flux through the ice must equal the atmospheric flux to the ice surface, as required by the energy balance (Drucker et al., 2003; Martin et al., 2004a). As in Martin et al. (2004a) and Drucker et al. (2003), we calculated the conductive heat flux, F_ice (W m^–2), through the ice assuming a linear temperature profile, with its lower boundary assumed to be at the seawater freezing point and its upper boundary an unknown ice surface temperature; as such, the conductive heat flux through the ice is given as,

where k_i = 2.03 W m^–1 K^–1 is the thermal conductivity of sea ice (Drucker et al., 2003), T_sea is the seawater temperature, assumed constant at –1.8°C, h_T is the derived ice thickness from Equation 3, and T_ice is the unknown ice surface temperature (Drucker et al., 2003).

Following Cavalieri and Martin (1994), we defined the total ice-to-atmosphere heat flux F_net just above the ice surface,

where F_T is the upward component of the turbulent heat flux, F_L is the longwave upward component, F_B is the longwave back-radiation component, and F_S is the shortwave downward component. F_T is defined as follows:

where ρ_air is the density of air (1.3 kg m^–3), C is the heat transfer coefficient (3.0 × 10^–3; Shokr and Sinha, 2015), c_p is the specific heat of air at constant pressure (1,004 J deg^–1 kg^–1), W is the 10-m air speed, and T_air is the 2-m air temperature. These latter 2 variables were obtained from the ERA5 reanalysis dataset (C3S, 2017) at 12.5-km spatial resolution. F_L is defined as follows:

where ε_ice is the emissivity of the ice surface (0.97; Shokr and Sinha, 2015) and σ is the Stefan–Boltzmann constant (5.6697 × 10^–8 W m^–2 K^–4). Finally, F_B and F_S were also obtained from the ERA5 reanalysis dataset, as “Surface thermal radiation downwards” and “Surface solar radiation downwards”, respectively.

To find the unknown ice surface temperature in Equation 4, we set Equations 4 and 5 equal to each other, then rearranged them so that all terms were on one side of the equation and set equal to 0. We then used an equation solver in Python to find the unknown ice surface temperature and then calculated the ice-to-atmosphere heat flux using Equation 5. A comparison of the derived ice surface temperature and air temperature shows the anticipated close relationship and dependency on ice thickness (Figure 3).

With heat flux calculated, we then calculated the total daily heat loss (H_L) in Joules from each grid cell defined as part of the polynya:

where 86,400 is the number of seconds in the day and A is the grid cell area (km²). The rate of freezing, dh/dt, in meters per second is given by:

where ρ_ice is the density of sea ice (0.92 × 10³ kg m^–3; Timco and Weeks, 2010) and L is the latent heat of fusion (3.34 × 10⁵ J kg^–1). A step further, the daily volume (km³) of ice production over the grid cell area was calculated as follows:

To demonstrate how this algorithm identifies a polynya and quantifies heat flux, we present an example of a large opening of the polynya on 5 February 2005 (Figure 4a). Figure 4a is a true color MODIS image that shows an area of new gray ice that had formed between coastal landfast ice and pack ice in central Hudson Bay, which appear white due to their thickness and potential snow cover. The fields of estimated ice thickness were overlaid in Figure 4b and show good agreement between the edge of thin ice within the polynya and the extent of new gray ice. The heat flux through the ice is presented in Figure 4c and highlights the gradient in heat flux, and subsequently heat loss and ice production, across the width of the polynya.

In Section 4, we examine the relationship between seasonal (December–April) average atmospheric conditions over northwestern Hudson Bay and annual cumulative ice production within the polynya. Daily gridded fields of sea level pressure, 10-m wind speeds, and 2-m air temperatures were downloaded from ERA5 (C3S, 2017) and used to calculate the seasonal averages. For a comparison of the relationship between northwesterly wind speed and ice production, the wind vectors contained within the study area (Figure 1) were projected at 49.16° such that they were perpendicular to the southeastern edge of the subregion (Figure 1); in this way, the offshoreward wind speed could be quantified and compared to ice production within the polynya.

Recent work has found improvement in ERA5 relative to its predecessor ERA-Interim (e.g., Betts et al., 2019; Graham et al., 2019); however, reanalysis data represent an assimilation of many datasets that are often sparse over the Arctic region due to its harsh climate. Thus, we compare the ERA5 wind speed data to in situ observations of wind speed at the Environment and Climate Change Canada meteorological station in Arviat (Figure 1). This station was selected because it had a complete record of hourly observations of wind speed and direction for the 16-year time period used in this study.

The primary objective of this study was to characterize the Kivalliq Polynya in terms of its spatiotemporal variability and ice production over 16 winters. The daily time series of cumulative ice production for each of the 16 years (Figure 5) reveals how variable and episodic ice production is within the polynya on an interannual basis. The 3 years with maximum ice production were 2004–2005, 2010–2011, and 2014–2015, and the 3 years with minimum ice production were 2003–2004, 2005–2006, and 2012–2013. The 16-year average of cumulative ice production was 182.2 km³. Generally, prolonged plateaus in ice production corresponded to periods when the polynya did not form and pronounced jumps corresponded with prolonged, large polynya events.

In the 16-year record, there were few days when the Kivalliq Polynya was not present in northwestern Hudson Bay, though its size was highly variable. Years with low ice production (2003–2004, 2005–2006, and 2012–2013) were characterized by relatively small openings and frequent closings (Figure 6). Conversely, years of high ice production (Figure 7) were characterized by persistent mid-sized polynya openings (e.g., 2004–2005 and 2014–2015) or the polynya remaining open for an extended period of time (e.g., 2010–2011). Indeed, the polynya closed only twice during this season: on March 16 and 17 and on April 8. Finally, 3 years were considered notable (Figure 8). The 2006–2007 and 2009–2010 winters had pronounced polynya events but were years with close to or below average cumulative ice production. Additionally, the 2013–2014 season showed relatively small but frequent openings similar in size to the low ice production years but was a year with above average cumulative ice production (Figure 5).

In general, there was a strong correlation between daily polynya area and ice production (R² = 0.660, P < 0.001), with coincident peaks and valleys in the two variables. This correlation was less pronounced at the beginning and end of the season during initial freeze-up and breakup of the ice cover in the area (i.e., in early December and late April). For example, in April 2006, there was a large increase in polynya area but a relatively smaller increase in ice production (Figure 6b), which was due to the seasonal increase in air temperatures that limited ice production within the existing polynya. In contrast, between freeze-up and breakup, when ice production in the region can be attributed to the opening of the polynya, there were both high wind speeds and low air temperatures that yielded rapid ice growth and subsequently a high peak in ice production. Examples of this situation can be seen in mid-February 2006 (Figure 6b), in early February 2005 (Figure 7a), and in mid-February 2011 (Figure 7b).

In contrast, there were occurrences throughout the 16-year record where a peak in polynya area did not lead to a similar peak in ice production, for example, in late January 2004 (Figure 6a), in early January 2007 (Figure 8a), and particularly in March 2010 (Figure 8b), which marked the beginning of a large, persistent polynya event and a relatively smaller increase in ice production. All of these events were marked by a peak in air temperature, suggesting a strong relationship between polynya area and air temperature. Over the 16-year record, there was a statistically significant correlation between daily air temperature and polynya area (R² = 0.140, P < 0.001) but not between daily air temperature and ice production (R² = 0.000, P = 0.375). At the seasonal scale, there was no statistically significant relationship between seasonally averaged air temperature and cumulative ice production (R² = 0.008, P = 0.734) or between seasonally averaged air temperature and annual median polynya area (R² = 0.005, P = 0.786).

Although the magnitude varies, the shape and location of the polynya remains similar interannually, with the highest number of polynya days concentrated in a long, narrow band along the western shore of Hudson Bay and declining with distance offshore into areas where only large polynya events extend (e.g., 2009–2010 and 2010–2011; Figure 9). Additionally, the polynya most commonly occurs in the central and southern end of the study region, south of Rankin Inlet, with limited activity in the northern portion of the study region near Roes Welcome Sound (Figure 9). In particular, 2009–2010 (a year with below average ice production) shows many polynya days concentrated along the coastline and in the southwest portion of the study area compared to the northeast portion (Figure 9).

Similar to the polynya occurrence maps in Figure 9, the ice production maps in Figure 10 show an offshore gradient with the greatest area of ice production focused along the landfast ice edge where frazil and thinner ice is typically present and, hence, the largest ocean–air heat fluxes occur. In particular, during high ice production years (2004–2005, 2010–2011, and 2014–2015), ice production is focused in this coastal band and extends throughout the study area. Conversely, during low ice production years (2003–2004, 2005–2006, and 2012–2013), the polynya does not extend throughout the full region and ice production in the coastal band is constrained. The greatest magnitude of ice production occurs south of Rankin Inlet with limited ice production in the northern end of the study area.

Because the Kivalliq Polynya is a latent heat polynya, variability in size and overall ice production between years can be linked to interannual variability in atmospheric forcing over the study area. Annual plots of the seasonal fields of sea-level pressure, surface winds, and air temperature reveal that years of high ice production (e.g., 2004–2005, 2010–2011, and 2014–2015) correspond to elevated sea-level pressure gradients and therefore stronger northwesterly winds over the study area, while years of low ice production (e.g., 2003–2004, 2005–2006, and 2012–2013) correspond to weaker sea-level pressure gradients and reduced northwesterly winds (Figure 11).

Statistically, we found that over the 16-year record, there is a significant relationship between daily average wind speed and daily polynya area with a one-day lag (R² = 0.078, P < 0.001) and a significant relationship between daily average wind speed and daily ice production (R² = 0.306, P < 0.001). These relationships imply that the instantaneous wind forcing does not strongly affect the size of the polynya but rather the cumulative wind effect appears to play a more important role in polynya area. There was a significant positive relationship (R² = 0.473, P = 0.003) between the seasonal average offshore wind speeds and cumulative ice production within the polynya that describes 47.3% of the variation in ice production (Figure 12). The slope of this linear regression indicates that ice production increases by 46.0 km³ for every 1.0 m s^–1 increase in seasonally averaged offshore wind speed (Figure 12). These observations are corroborated by in situ observations of surface winds from Arviat, Nunavut, that show an average seasonal northwesterly wind during winter (red arrow; Figure 11) and a significant positive correlation with the ERA5 (zonal: R² = 0.740, P < 0.001; meridional: R² = 0.726, P < 0.001). Regardless of the interannual variability in wind speed magnitude, higher wind speeds occurred over the part of the study area south of Rankin Inlet than north of it (Figure 11), which was likely a contributing factor to the higher magnitude of ice production south of Rankin Inlet seen in Figure 10 for each year.

As the Kivalliq Polynya contributes significantly to the Hudson Bay ice mass balance, examining the interannual variability of ice production in the polynya and how this variability may be explained by variation in wind forcing and heat flux is important. Cheng et al. (2019) conducted a quantitative analysis of this type on the Ross Ice Shelf Polynya by examining the correlation between ice production and polynya area, wind speed, and heat flux components. Their findings indicated that ice production is strongly correlated to polynya area. Likewise, from our results, polynya frequency (i.e., those grid cells with ice thickness 10 cm or less; Figure 9) directly relates to ice production (Figure 10), with a consistent narrow band of frequent polynya events and high ice production along the northwest coast of Hudson Bay and, in most years, a strong shift toward cells with fewer open water days/lower ice production away from the coastline.

As noted in Section 3, there was a significant positive correlation between seasonal average wind speed and cumulative ice production, such that a 1.0 m s^–1 increase in wind speed results in a 46.0 km³ increase in ice production (Figure 11). High ice production in 2004–2005, 2013–2014, and 2014–2015 corresponded to a strong sea-level pressure gradient over northwestern Hudson Bay and subsequently high northwesterly wind speeds; similarly, low ice production in 2003–2004, 2005–2006, and 2012–2013 corresponded to a weak sea-level pressure gradient and low wind speeds (Figure 11). Beyond these high and low ice production years, the remaining years also had weak sea-level pressure gradients and low wind speeds. For example, 2009–2010 had slightly below average cumulative ice production but a large, extended polynya event that began at the end of February; this event could be expected to be the result of high wind speeds that are not reflected at the seasonal scale. However, plots of daily ice production (Figure 8) make apparent that warm air temperatures rather than wind speed played a greater role in this event. Conversely, high ice production in 2010–2011 was far greater than what would be expected based on the seasonal average wind speed, but as depicted in Figure 7, steady northwesterly winds persisted throughout winter 2011 and facilitated the extended opening of the polynya that led to increased ice production.

Both Morales Maqueda et al. (2004) and Smith and Barber (2007) suggested that because of their transitory nature, polynyas are useful indicators of large-scale climate change. However, the effects of climate change can have confounding impacts on a polynya. Increasing air temperatures not only affect heat flux through the ice but have also shortened the duration of the ice cover throughout Hudson Bay (Steiner et al., 2013; Andrews et al., 2017) and thereby shortened the period between freeze-up and break-up when the polynya can exist. Indeed, between 2020 and 2070, northwestern Hudson Bay is projected to warm between 0.06°C and 0.13°C on average per decade for RCP 4.5 and between 0.12°C and 0.27°C on average per decade for RCP 8.5 (P Myers, personal communication). However, ice production in the polynya is driven predominately by prolonged or episodic openings that are driven by offshore winds; hence, the effect of climate change on the prevailing northwesterly winds over Hudson Bay will exert a much greater influence on the occurrence and size of the polynya than will increased air temperature. In particular, our analysis showed that ice production responds to small changes in offshore wind speeds, suggesting that small changes in the magnitude of the northwesterly winds will have a considerable impact on the polynya and ice mass balance of Hudson Bay. These potential changes would impact primary production within northwestern Hudson Bay and the marine mammals that rely on the food web of this area (e.g., the western Hudson Bay seal and polar bear population; Smith and Barber, 2007). As such, the analysis we have conducted of the spatiotemporal variability of the Kivalliq Polynya provides some guidance on how the polynya may respond to climate change and may be useful as an indicator of climate change in the Hudson Bay region.

The largest inland sea in the world, Hudson Bay is isolated from the open water circulation of the Arctic and Atlantic Oceans; thus, variations in Hudson Bay’s seasonal ice cover are primarily the result of atmospheric forcing. Prevailing northwesterly winds advect the pack ice across the Bay, creating thicker deformed ice in eastern Hudson Bay and maintaining the Kivalliq Polynya in northwestern Hudson Bay. Using 16 years of daily passive microwave data and a thin ice algorithm, we have provided the first detailed analysis of the spatiotemporal variability of this regionally important polynya during winter. On average, 182 km³ of sea ice are produced within the polynya annually, though the interannual variability is considerable. Large, persistent polynyas through the winters of 2004–2005, 2010–2011, and 2014–2015 led to anomalously high ice production, whereas the polynya remained relatively small through the winters of 2003–2004, 2005–2006, and 2012–2013, leading to anomalously low ice production during these years. Spatially, ice production is focused in the southern end of the study region, south of Rankin Inlet, and concentrated in a narrow band along the landfast ice edge of northwestern Hudson Bay before declining with distance offshore.

We found a significant relationship between the average seasonal offshore wind speeds and cumulative ice production that corresponds to an increase of 46.0 km³ in ice production for every additional increase of 1.0 m s^–1 in offshore wind speeds. Broadly, years with weak offshore winds are characterized by smaller polynyas and less ice production (e.g., 2003–2004, 2005–2006, and 2012–2013), whereas years with pronounced offshore winds have large polynyas and greater ice production (2004–2005, 2010–2011, and 2014–2015). However, some years are characterized by large episodic polynya events, which are driven by periods of pronounced offshore winds, but become outliers in the seasonal relationship between offshore winds and ice production (e.g., 2009–2010). At a daily timescale, these outliers can be explained by a stronger relationship between ice production and air temperature rather than ice production and wind speed.

The datasets used in this analysis are publicly available and listed in the Reference section: C3S (2017), Cavalieri et al. (2014a, 2014b), Markus et al. (2018), Meier et al. (2018), and Meier et al. (2019).

The authors would like to thank Kevin Sydor and Karen Wong from Manitoba Hydro for their input during manuscript revision. Thank you to the National Snow and Ice Data Center, which provided the SSM/I, AMSR-E, and AMSR-2 brightness temperature data, and to the Copernicus Climate Change Service for the provision of the ECMWF ERA5 reanalysis dataset. The authors would also like to thank Geoffrey Gunn, who created the shapefile in Figure 1.

This work is a part of the BaySys program examining the Hudson Bay system, funding for which is provided by Manitoba Hydro. D. Babb, J. Ehn and D.G. Barber are supported by the Natural Sciences and Engineering Research Council of Canada. D. Babb is additionally supported by the Canadian Meteorological and Oceanographic Society. Financial support was also provided by the Canada Excellence Research Chair and the Canada Research Chair programs. This is a contribution to the ArcticNet and Arctic Science Partnership networks.

The authors declare no conflict of interest regarding the publication of this manuscript.

Contributed to conception and design: JB, DB, WC, DGB.

Contributed to acquisition of data: JB, WC.

Contributed to analysis and interpretation of data: JB, DB, WC.

Drafted and/or revised the article: JB, DB, WC, SK, JE, JH, DGB.

Approved the submitted version for submission: All authors.

How to cite this article: Bruneau, J, Babb, D, Chan, W, Kirillov, S, Ehn, J, Hanesiak, J, Barber, DG. 2021. The ice factory of Hudson Bay: Spatiotemporal variability of the Kivalliq Polynya. Elementa: Science of the Anthropocene 9(1). DOI: https://doi.org/10.1525/elementa.2020.00168

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Stephen F. Ackley, Department of Geological Sciences, University of Texas at San Antonio, TX, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: BaySys