IcePAC – a probabilistic tool to study sea ice spatio-temporal dynamics: application to the Hudson Bay area

A reliable knowledge and assessment of the sea ice conditions and their evolution in time is a priority for numerous decision makers in the domains of coastal and offshore management and engineering as well as in commercial navigation. As of today, countless research projects aimed at both modelling and mapping past, actual and future sea ice conditions were completed using sea ice numerical models, statistical models, educated guesses or remote sensing imagery. From this research, reliable information helping to understand sea ice evolution in space and in time is available to stakeholders. However, no research has, until present, assessed the evolution of sea ice cover with a frequency modelling approach, by identifying the underlying theoretical distribution describing the sea ice behaviour at a given point in space and time. This project suggests the development of a probabilistic tool, named IcePAC, based on frequency modelling of historical 1978–2015 passive microwave sea ice concentrations maps from the EUMETSAT OSI-409 product, to study the sea ice spatio-temporal behaviour in the waters of the Hudson Bay system in northeast Canada. Grid-cell-scale models are based on the generalized beta distribution and generated at a weekly temporal resolution. Results showed coherence with the Canadian Ice Service 1981–2010 Sea Ice Climatic Atlas average freeze-up and melt-out dates for numerous coastal communities in the study area and showed that it is possible to evaluate a range of plausible events, such as the shortest and longest probable ice-free season duration, for any given location in the simulation domain. Results obtained in this project pave the way towards various analyses on sea ice concentration spatio-temporal distribution patterns that would gain in terms of information content and value by relying on the kind of probabilistic information and simulation data available from the IcePAC tool.

1 Introduction

Numerous scientific projects have recognized the link between climate change and changes in the spatio-temporal sea ice distribution (Andrews et al., 2017; Cavalieri and Parkinson, 2012; Comiso, 2011, 2002; Comiso et al., 2008; Gloersen et al., 1999; Johannessen et al., 2004; Rothrock et al., 1999; Stocker, 2014; Stroeve et al., 2007; Stroeve et al., 2014; Stroeve et al., 2012; Wang and Overland, 2009). The climate change in the Arctic (60–90^∘ N) is expected to be, on average, 1.9 times greater than the global mean warming (Winton, 2006), and this “Arctic amplification” is expected to strengthen in coming decades (Serreze and Barry, 2011).

In this changing environment, an adequate and efficient monitoring of sea ice is of key importance in better understanding the climate and its impacts on marine and coastal areas (Barnhart et al., 2014; Bintanja and Selten, 2014; Davies et al., 2014; Holland et al., 2006; Kowal et al., 2017; Manabe and Stouffer, 1995; Overeem et al., 2011; Peterson et al., 2002; Rahmstorf, 1995; Rahmstorf and Ganopolski, 1999; Vermaire et al., 2013), on the security of economical and logistical activities in northern communities (Aksenov et al., 2017; Andrews et al., 2017; Ho, 2010; Lasserre and Pelletier, 2011; Liu and Kronbak, 2010), on arctic marine fauna protection (Bhatt et al., 2010; Castro de la Guardia et al., 2013; Darnis et al., 2012; Laidre et al., 2015; Post et al., 2013; Wassmann et al., 2011), and on the traditional way of life of Inuit communities (Durkalec et al., 2015; Laidler et al., 2010).

To understand and appreciate the role of sea ice cover regarding climate, marine and coastal environment management, fauna protection, and the cultural traditions of northern communities, access to informative and reliable sea ice spatio-temporal distribution information is fundamental. Engineers, stakeholders, Inuit and northern populations, navigators, and scientists must be able to quantify hazards related to the sea ice cover in order to efficiently evaluate, anticipate, and minimize the risks of usage, building, and exploitation in marine and coastal areas. Given the increase in activity like tourist cruises, shipping, and mining observed in the Arctic and the north (Dawson et al., 2018; Lasserre and Pelletier, 2011; Pizzolato et al., 2016), one can expect the demand in information to also increase. For example, engineers could make use of probabilistic data to assess the potential duration of sea ice presence for an infrastructure they are planning to build; mariners could use the data to estimate the best departure date from their home port to reach their final destination according to a certain sea ice concentration probability threshold; fauna specialists could use the data to estimate the risk encountered by species dependent on sea ice cover for their fitness, such as polar bears and seals; and finally, Inuit communities could use the tool to evaluate if their planned travel routes are risky for a given period of the year given the known history of the sea ice spatio-temporal behaviour.

Despite a large number of Earth observation datasets on sea ice cover, only a few provide both high temporal and spatial resolution. National ice services, such as the Canadian Ice Service (CIS), provide users with sea ice conditions' climatology that is a reliable source of descriptive statistics on the sea ice spatio-temporal behaviour, such as the average freeze-up or maximum extent date. The CIS also provides daily and weekly regional ice condition maps that inform users of the observed concentrations and the ice development stage reached by sea ice, as well as detailed sea ice condition reports, for all regions of the Canadian Arctic. These maps and reports are created by experienced and skilled professional sea ice analysts who make use of diverse sources of sea ice information such as radar, optical microwave imagery, and passive microwave imagery in combination with in situ observations to prepare their analyses (Iacozza, 2000). However, even if the CIS data and other national ice service products do provide probabilistic information, these datasets do not carry information on the nature of the underlying statistical distributions of sea ice parameters, such as sea ice concentration (SIC %), at any given point.

To build a probabilistic model of sea ice concentrations, historical information is needed and must meet specific needs such as long-term availability, reliability, large spatial coverage and high temporal frequency. As visible imagery is largely affected by cloud coverage, especially prominent in the Arctic, it is not a reliable source of information with which we can build our model. In spite of its independence of atmospheric conditions, synthetic aperture radar (SAR) imagery does not provide sufficient spatial coverage for our purpose. Therefore, passive microwave observation turned out to be a compromise as this dataset, even if its resolution is coarse, provides daily data for the entire Arctic and historical data are available for it from 1978 onwards.

By exploring an innovative probabilistic sea ice concentration modelling avenue, this study proposes a tool, named IcePAC, to characterize the underlying statistical distributions of the SIC % at any point in the Hudson Bay system (Saucier et al., 2004) based on historical passive microwave remote sensing data from 1978 to 2015. These data are than used to analyze the spatio-temporal behaviour of SIC % in the Hudson Bay area, with a probabilistic perspective and compared to the CIS climatology.

Figure 1 The Hudson Bay system with communities, model validation points and polynyas.

2 The Hudson Bay system

The study area is the Hudson Bay system (HBS), consisting of the Hudson Bay, Hudson Strait, James Bay, and Foxe Basin (Fig. 1). The HBS is surrounded by the three Canadian provinces of Quebec, Ontario, and Manitoba and the territory of Nunavut. It is the largest inland sea on Earth, with a total area of 1 300 000 km² (Etkin, 1991; Gagnon and Gough, 2005a; Martini, 1986), and is located in both subarctic and Arctic regions. An estimated 20 % of the flux of freshwaters to the Arctic Ocean are thought to come from rivers flowing into the HBS, which represents 900 km³ yr⁻¹ (Déry and Wood, 2005, 2004). It is connected to the Labrador Sea via the Hudson Strait while waters from the Arctic Ocean flow through the Fury and Hecla Strait to the Foxe Basin (Prinsenberg, 1986), and it is characterized by shallow depths of less than 100 m in the Foxe Basin, of a maximum of 125 m in the Hudson Bay, and of more than 200 m in the Hudson Strait (Jones and Anderson, 1994). It has cyclonic currents generated mostly by winds, with a maximum intensity in November (Saucier et al., 2004).

A large amount of research has been done to document the average sea ice behaviour in the HBS (CIS, 2013; Gagnon and Gough, 2005a, b; Hochheim and Barber, 2010, 2014; Kowal et al., 2017; Maxwell, 1986), which goes through a complete freeze–thaw cycle every year. The sea ice cover in the HBS is primarily constituted of first-year ice, with the exception of traces of multi-year ice drifting in the Foxe Basin (CIS, 2013; Etkin and Ramseier, 1993). The sea ice cover initially forms in the northwestern part of the HBS near Southampton Island in late November and progresses towards the southeastern part of the HBS (Hochheim and Barber, 2014; Maxwell, 1986) to finally cover most of the HBS in late December. Sea ice maximum extent in the HBS is usually achieved in April (Gagnon and Gough, 2005b), after which the melt begins in May along the northwestern shoreline of the HBS. The melt progresses from the shores toward the centre of the Hudson Bay, which usually results in an agglomerate of sea ice in the south central part of the Hudson Bay in late July (CIS, 2013). In summary, the HBS is, on average, frozen in late December and free of ice in mid-August (Mysak et al., 1996; Wang et al., 1994).

Many authors have studied the trends in sea ice cover extent for the Hudson Bay area (Galbraith and Larouche, 2011; Hochheim and Barber, 2014; Tivy et al., 2011). Among them, Tivy et al. (2011) arrived at the conclusion that the Hudson Bay area was affected by some of the strongest downward trends regarding the ice season duration in the entire circumpolar Arctic. Confirming the results of Tivy et al. (2011), Hochheim and Barber (2014) measured trends of the open water season duration comparing a 1996–2010 climatology with a 1980–1995 climatology based on a modified Comiso SIC % dataset. Their results showed lengthening of ice-free seasons in Foxe Basin, Hudson Strait, and Hudson Bay, of 3.5, 4.9, and 3.1 weeks respectively. It is worth noting that, since the HBS is mostly covered by first-year ice, the natural variability of the sea ice conditions is considerable since it is mostly driven by warming temperatures, but also by changes in atmospheric circulation (Mudryk et al., 2018).

Ice thickness during winter in the HBS ranges from 1 to >2.5 m according to numerical modelling studies, though, as reported by Landy et al. (2017), these studies do not agree on the spatial distribution of sea ice. Gough et al. (2004) identified an east–west asymmetry in long-term trends of sea ice thicknesses in the HBS using drill-hole measurements acquired between 1960 and 2000. These trends show a tendency of thickening on the western side (+0.1 to 1.5 cm yr⁻¹) of the Hudson Bay, while the eastern side shows, conversely, a trend towards a thinning ice pack (−0.5 to 0.8 cm yr⁻¹).

Polynyas are also present in the HBS, such as the northwestern Hudson Bay polynya between the western coast of the Hudson Bay and Southampton Island that forms occasionally throughout the winter and spring (Gough et al., 2004), the Fury and Hecla Strait and Hall Beach polynyas (Barber and Massom, 2007), located in northwestern parts of the Foxe Basin, and the Cape Dorset polynya, which is a shore lead polynya (Stirling, 1980).

3 Data and methods

Sea ice extent has displayed an important decline in the last decades (Cavalieri and Parkinson, 2012), as it can be observed with remote sensing data, such as passive microwave data, which have been acquired since 1978. Another important source of information on sea ice cover is model predictions which come from deterministic models (Hunke et al., 2017; Rousset et al., 2015; Weaver et al., 2001), based on dynamical and thermodynamic equations evolving in synergy inside a modelling framework, or from statistical models, based on statistical tools such as simple and multiple regression analysis (Ahn et al., 2014; Drobot, 2007; Pavlova et al., 2014) to explain an expected sea ice parameter value (e.g. sea ice extent, sea ice area, sea ice concentration, sea ice thickness).

Another statistical approach, focusing on the estimation of the probabilities of occurrence of specific sea-ice-related events, has been used in recent research (Dirkson, 2017; Rajak et al., 2015). It is achieved either by using the simple count method (e.g. an event occurred four times in the last 10 years, which corresponds to a 40 % probability of occurrence) or by using the frequency modelling method, which consists of adjusting a theoretical distribution to a series of observations, consequently defining the plausible sea ice events for the entire range of probabilities (i.e. p=0 to 100 %). In this research, the frequency modelling method is used for a series of passive microwave historical SIC % remote sensing data to adjust distributions to a total of 20 738 grid cells or sites within the HBS (i.e. the spatial dimension) for each of the 52 weeks of the year (i.e. the temporal dimension), resulting in a total of 1 078 376 distribution fits.

The datasets used and protocols followed in the IcePAC tool to model SIC % distributions at every grid cell in the HBS are described in the following sections.

3.1 Sea ice concentration dataset

Sea ice concentration is defined as the proportion of sea ice covering a predefined area, expressed as a percentage. In remote sensing, this predefined area is represented by a grid cell. The choice of a SIC dataset for frequency modelling in our study is highly influenced by the extent of the HBS (1 300 000 km²) and has been made to ensure uniformity and continuity of the series used for analysis.

Multiple SIC datasets are generated using either visible, thermal, SAR, or passive microwave remote sensing. As the objective of the IcePAC tool is to provide the capacity to evaluate the spatio-temporal evolution of the SIC, the source of data needed to ensure a complete coverage of the HBS (i.e. the spatial dimension) and to ensure continuity in the series (i.e. the temporal dimension). A passive microwave dataset was chosen as it meets these two needs.

The global reprocessed sea ice concentration dataset, OSI-409 (Eastwood et al., 2015; Tonboe et al., 2016), was selected as it enables the reconstitution of SIC series for more than 30 years with a 12.5 km grid size and is processed with a unique hybrid SIC algorithm. The hybrid algorithm only uses the information taken from the “Bootstrap” algorithm (Comiso, 1995) when SIC <70 %, linearly weights the SIC estimated by the Bootstrap and Bristol (Smith, 1996) algorithms when 70 % < SIC > 90 %, and only uses the information from the Bristol algorithm when SIC > 90 %. Another passive-microwave-based SIC % dataset is used in this study as a comparison dataset, the OSI-430, and its only difference from OSI-409 is that it uses SSM/I data obtained from NOAA instead of recalibrated SSM/I data from RSS (Remote Sensing Systems). The difference in the resulting SIC % product is, according to Eastwood et al. (2015), expected to be minimal.

In OSI-409, the passive microwave channels used for ice concentration mapping have footprint sizes ranging from 56 km for the 19 GHz channels to 33 km for the 37 GHz channels. As SIC % values are represented on 12.5 km grid cells in the OSI-409 and 430 products, inputs of different resolutions are combined using a gridding procedure that loads all passive microwave observations within the period of a day for a 12.5 km grid cell (centred on 12:00 UTC) and averages them using a weighting value (dependent on the distance between the observation and the centre of the grid cell) and an influence radius (dependent on the passive microwave channel resolution). A detailed explanation of the method is provided in Eastwood et al. (2015).

It is important to note that passive microwave SIC % datasets are known to be affected by diverse error factors such as a land spill-over effect along coasts that triggers false higher SIC % estimation if not taken into account. To mitigate the errors caused by this phenomenon, a coastal correction is applied to the data with a method inspired from Cavalieri et al. (1999). This method first calculates monthly average SIC % for all months and finds the minimum ice concentration from these averages. This minimum is then used to correct the ice concentration values in the coastal zone if adjacent non-coastal grid points are ice-free. Also, climatological maximum extent masking is done to mask out erroneous ice outside areas where sea ice is ever likely to occur using a sea ice extent monthly climatology from the NSIDC (National Snow and Ice Data Center) (NSIDC, 2013).

Figure 2 The process of frequency modelling from series building to model fit. (a) The SIC % data stacking and series building for every pixel in the simulation domain. (b) An example of a unique site complete series plotted. (c) An example of model fit (beta) to the density distribution of the SIC % values in the series with the α and β model parameters.

Another error factor with passive microwave data is data gaps which can occur either in the form of missing scan lines, missing orbits, or polar observation holes. As reported by Eastwood et al. (2015), these gaps, when small, can be corrected through simple interpolation. However, when facing large gaps, a blurring effect appears in the interpolated area. To correct this effect, an interpolation approach using the information from the past and following days is used in the OSI products. And last but not least, the effect of ponds appearing during the melt is that the resulting maps tend to underestimate SIC % during summer since there is confusion between open water areas and melt ponds on top of sea ice. Finally, another underestimation of the SIC % results from thinner ice types which do not act as a radiometric insulator for the passive microwave frequencies around 19 and 37 GHz that are the base of the OSI-409 and OSI-430 datasets (Eastwood et al., 2015).

These different error sources in the process of estimating SIC % using passive microwave imagery are known to have an impact on the reliability of the data, especially during the freeze-up and the melt periods, as brought forward by Agnew and Howell (2003), who noticed that the underestimation of ice extent in the Hudson Bay during summer, when compared to CIS maps, could go up to 43.5±27.9 % (while considering the CIS data as ground truth). Since the HBS is an area where new ice forms every year, this systematic underestimation when using passive microwave data must be kept in mind.

The data have been clipped to the HBS extent using the Natural Earth (NaturalEarth, 2014) vector dataset with an estimated spatial resolution of 500 m (Wessel and Smith, 1996), well beyond the resolution of the OSI-409 and OSI-430 datasets. As the coasts of the HBS are highly dynamic, other datasets could have been used such as the CanVec product from Natural Resources Canada, which is updated regularly. However, considering the 12.5 km grid size, the Natural Earth product was chosen.

3.2 Frequency modelling

The IcePAC tool uses frequency modelling to describe the underlying SIC % probability distribution at a given site with a simplified model fitted on historical SIC data. The first step in this approach is to build the time series of historical data, then to ensure their quality using preliminary tests, and finally to identify and fit the model on the data (Fig. 2).

3.2.1 Building the sea ice concentration (SIC) series

The SIC series are built to represent the SIC state for a specific week, for all years between 1978 and 2015. It is for these series that we adjust a theoretical distribution to estimate the probabilities of SIC-related events. First, the daily OSI-409 data have to be averaged every 7 days to create weekly datasets. This operation is made following a 365-day “no-leap” calendar convention (i.e. every year has 365 days), separated in 52 weeks (31 December is included in week 52). Second, the data for each week number are stacked in chronological order, from 1978 to 2015.

Figure 3 The trend removal process applied to a SIC % time series for a coastal grid cell extracted from OSI-409 data for week 50.

3.2.2 Preliminary tests

Series have to go through a set of preliminary tests to assess their stationarity, homogeneity, and independence, assuring they are suitable for frequency modelling. The tests used in IcePAC are the Mann–Kendall test for stationarity (Mann, 1945), the Wald–Wolfowitz test for independence (Wald and Wolfowitz, 1940), and the Wilcoxon test for homogeneity (Wilcoxon, 1945).

Series are considered independent if the subsequent unique SIC observations have no dependence on one another, they are considered homogeneous if they are reputed to be from the same distribution, and they are stationary if they are not affected by a trend. In the case of a detected non-stationarity, the trend is modelled and subtracted from the series (Cave and Pearson, 1914).

All the time series, either original or detrended in regard to the Mann–Kendall test, satisfied the preliminary tests for a significance level of α=0.05.

3.2.3 Trend estimation and removal

The estimation of a trend on a percentage data series has the particularity that the trend must, in order to be coherent with the physics of the studied phenomenon, be bounded in a [0,1] domain (Baum, 2008). In other words, we must ensure we do not measure a trend that generates SIC values larger than 100 % or smaller than 0 %. To guarantee the respect of this criterion, a generalized linear model with a logit link function has been used to estimate the trend. The logit link function linearizes the SIC values using the logit transformation (Eq. 1) and it is with these transformed values that a linear regression of the form $\mathit{\alpha }x+\mathit{\beta }+\mathit{\epsilon }$ is measured.

$$\begin{array}{}\text{(1)}& \text{logit}\left(\text{SIC}\right)=\ln (\text{SIC}/\mathrm{1}-\text{SIC}),\end{array}$$

where SIC is defined ]0,1[.

For the trend to be removed from the series, an inverse transformation (Eq. 2) must be applied to the estimated trend to turn its logit values into SIC values.

$$\begin{array}{}\text{(2)}& \text{SIC}=\exp \left({\displaystyle \frac{\exp \left(\text{logit}\left(\text{SIC}\right)\right)}{\left(\mathrm{1}+\exp \left(\text{logit}\left(\text{SIC}\right)\right)\right)}}\right)\end{array}$$

The removal of the trend in the time series ensures that we model the natural variability of the sea ice concentration phenomenon, without any influence from the trend (Fig. 3). In IcePAC, the trend is modelled using the aforementioned method (GLM with logit link), then it is removed from the original SIC % data, and finally the residuals are used to adjust the distribution model.

Once adjusted, the distribution models are used for queries, in conjunction with the trend that is taken into account to generate the final result. For example, one could ask for the probability of a specific SIC %. In such a case, the trend for the prediction year is first removed from the SIC %, and the probability of the residual is then estimated from the distribution model. Inversely, one could ask for the SIC % for a given probability. In such a case, the probability is first used to get the corresponding residual value (representing the natural variability of sea ice) and the trend is reinjected afterward to obtain a realistic SIC % value. In both cases, the trend is taken into account when generating the final results.

It is important to note that frequency analysis is not a projective approach but a predictive approach. In other words, IcePAC is not to be used to get an outlook of long-range future SIC % conditions (projection) but to assess what is expected for the short-term conditions (prediction). As the trend will change in time with the addition of new SIC % data, the further we try to temporally expand our prediction, the more erroneous our probability estimates may be.

3.2.4 Distribution selection and fit

The selection of adequate candidate distributions to fit on the series is largely limited by the bounded nature of the data. The selected distribution has to be bounded to [0,1] and be available in a generalized form in order to adapt to detrended series which are bounded to [−1,1]. Their generalized forms are to be used with the position parameter a fixed at −1 and the scale parameter b fixed at 2, in coherence with the phenomenon. It must also present enough flexibility in shape to adapt to the different type of SIC series in the HBS domain. Two different distributions, the generalized beta distribution and Johnson's SB (system-bounded) distribution (Johnson, 1949), both bounded and displaying flexibility in shape, have been fitted on the series using the maximum likelihood estimator (MLE; NIST, 2013) and compared by measuring the root mean square error (RMSE) between observations and adjusted curves as well as the Akaike information criterion (Akaike, 1998) and the Bayesian information criterion (Schwarz, 1978).

The generalized beta distribution (Eq. 3) has four parameters which are the p (p>0) and q (q>0) shape parameters, the position parameter $a=-\mathrm{1}$, and the scale parameter b=2. In Eq. (3), B is the beta function. This distribution has been used before in climatology (Dirkson, 2017; Falls, 1974; Henderson-Sellers, 1978; Sulaiman et al., 1999), in seismology (Lallemant and Kiremidjian, 2015), to study air pollution (Nadarajah, 2008) and in hydrology (Chen and Singh, 2017; Yao, 1974).

$$\begin{array}{}\text{(3)}& f\left(x\right)={\displaystyle \frac{(x-a{)}^{p-\mathrm{1}}(b-x{)}^{q-\mathrm{1}}}{B(p,q)(b-a{)}^{p+q-\mathrm{1}}}},\end{array}$$

where $a\le x\le b$; $p,q>\mathrm{0}$.

Johnson's SB (Eq. 4) is, under its generalized form, a four-parameter distribution, for which the γ and δ (δ>0) are the shape parameters, the parameter $a=-\mathrm{1}$ is the position, and the parameter b=2 is the scale. This distribution has been used before in meteorology (Cugerone and Michele, 2015; Wakazuki, 2013), in forestry (Rennolls and Wang, 2005), and in hydrology (D'Adderio et al., 2016).

$$\begin{array}{}\text{(4)}& f\left(x\right)={\displaystyle \frac{\mathit{\delta }}{\sqrt{\mathrm{2}\mathit{\pi }}}}{{\displaystyle \frac{(b-a)}{(x-a)(b-x)}}}^{\left[-\frac{\mathrm{1}}{\mathrm{2}}{\left\{\mathit{\gamma }+\mathit{\delta }\ln \left(\frac{x-a}{b-x}\right)\right\}}^{\mathrm{2}}\right]},\end{array}$$

where $b,\mathit{\delta }>\mathrm{0}$.

Two approaches were tested for distribution selection. The first approach considered fitting the distributions to the series according to their Mann–Kendall test results (i.e. managing detected trends only). The second approach considered fitting the distributions to systematically detrended series (i.e. removing the trend from every time series in the simulation domain) in order to ensure spatial coherency in the IcePAC outputs.

The use of a systematic trend removal approach is justifiable by the fact that natural phenomena are considered by nature non-stationary (Lins, 2012; Rao et al., 2012) and by the relative shortness of the series which can have an effect on the conclusions of the Mann–Kendall test (Hirsch et al., 1982). Also, it is important to state that the frequency analysis approach is rarely used to generate spatialized results like it is in the IcePAC approach in which every time series (linked to a specific location or grid cell) is processed as an individual station.

The two trend removal approaches, tested on 958 randomly chosen series, yielded similar conclusions. In 71.4 % of cases, the beta distribution outperformed Johnson's SB distribution. For the remaining 28.6 % adjustments for which Johnson's SB distribution did perform better, both the information criterion and the RMSE showed a non-significant difference with the beta distribution. In light of these results and to preserve parsimony in IcePAC, it was decided to use only the beta distribution for all series, for which a systematic trend removal was applied.

As expected, using an approach which systematically corrects for trends in time series before the distribution fit improved the spatial coherency of the resulting probability maps generated by IcePAC (Fig. 4).

Figure 4 Comparison of IcePAC results for non-systematic (i.e. based on Mann–Kendall test result) trend removal (a) and systematic trend removal (b), on the resulting map for the probability of observing a SIC < 50 % on week 1.

3.2.5 Model queries

Queries with the IcePAC tool are possible via three important functions resulting from the distribution fit, the probability density function (PDF), the cumulative distribution function (CDF), and the percent point function (PPF). These functions are obtained using the fitted parameters from the beta distribution independently for each of the 20 738 grid cells in the IcePAC simulation domain.

The PDF is obtained by fitting a theoretical distribution on the frequency histogram of the SIC % observations. The selected distribution in IcePAC is the beta distribution for which the parameters of Eq. (3) were estimated using the MLE. The PDF shows how probability is distributed between SIC % values and how it evolves.

Derived from the PDF, the CDF (Eq. 5) gives the probability for a given range of SIC % values. It corresponds to the area under the PDF curve for a specific range of SIC values, usually from 0 % to SIC %_max. As an example of the CDF, one could query the probability of non-exceedance (p) for a sea ice concentration of 25 % (SIC_max) for week number 1 for the year to come.

$$\begin{array}{}\text{(5)}& F_x\left(\text{SIC}\mathit{\%}\right)=\underset{\mathrm{0}}{\overset{\text{SIC}{\mathit{\%}}_{\text{max}}}{\int }}{f}_x\left(\text{SIC}\mathit{\%}\right)\text{d}t\end{array}$$

The inverse function of the CDF, the PPF (Eq. 6), estimates the SIC % value for a given probability of non-exceedance. As an example of PPF, one could query the SIC %_max for a probability of non-exceedance of 55 % for week 1 for the year to come.

$$\begin{array}{}\text{(6)}& Q=F^{-\mathrm{1}}\left(p\right)\end{array}$$

Since the IcePAC fits are made on detrended series (e.g. residuals), the trend has to be taken into account when processing queries, meaning that it has to be either removed from the SIC %_max value if the CDF is used or added to the result of the PPF query in order to render a physically valid result. The query flow chart is presented in Fig. (5).

3.2.6 IcePAC versus observations in 2016

The assessment of the validity of IcePAC predictions was done by comparing IcePAC weekly output time series for the entire year with the OSI-430 product, a data source not used in IcePAC development but based on the same SIC % retrieval algorithm (Tonboe et al., 2016). The comparison was made between the 2015–2016 sea ice season SIC % values (not included in the input data for IcePAC) and the IcePAC SIC % for a non-exceedance probability of 90 % (P=0.9)

Eight different comparison sites were selected to represent different sea ice spatio-temporal behaviours (see Fig. 1). Four coastal sites, Akismi Island (CAI), Cape Dorset (CCD), Belcher Islands (CBI), and Hall Beach (CHB), were sampled to assess the behaviour of IcePAC predictions along the coastline at different latitudes. Also, four offshore sites, Frobisher Bay (OFB), Central Hudson Bay (OCHB), Churchill (OC), and Northern Ungava Bay (ONUB) were sampled to assess the behaviour of IcePAC predictions offshore at different latitudes and at critical navigation passage points.

Figure 5 IcePAC query flow chart.

Figure 6 displays SIC % prediction outputs from IcePAC, for a probability of non-exceedance of 90%. Each comparison site shows a modelled dynamic which is coherent with reality and with the OSI-430 2015–2016 observation data. As it could be expected, for all sites, it is during the freeze-up and melt periods that we can observe the largest differences between the 2015–2016 OSI-430 SIC % observations and the model outputs, compared to stable cover periods during which the range of probable SIC % values is narrowed down (i.e. the mean and the P=90 % lines are almost overlaid), and therefore the observations tend to rejoin with the model output. Some anomalies of various intensities can be detected in the series such as early 2015 melt-out events at Hall Beach and Churchill (▵ in Fig. 6), late 2016 freeze-up events at points Akismi Island and Churchill (∘ in Fig. 6), and early 2016 melt-out event at Hall Beach (◊ in Fig. 6), which all are in agreement with the anomaly maps (Fetterer et al., 2017) of the NSIDC's Arctic Sea Ice News and Analysis (https://nsidc.org/arcticseaicenews/, last access: 18 May 2018). The fact that the Cape Dorset site displays a SIC % for a probability of non-exceedance of 90 % of about 80 % at most is coherent with its shore lead polynya status. Similar conclusions can be made for the Hall Beach point for which the polynya tends to appear by late May, which in our case was earlier in 2015 and 2016, as confirmed by NSIDC. In Fig. 6, daily model outputs (green lines) are presented as a comparison with the weekly outputs (black dots) in order to justify why the IcePAC model outputs were generated at a weekly interval, to filter the effects of the daily variability of SIC % estimations in OSI-409, especially visible at coastal sites.

Figure 6 IcePAC weekly and daily P=90 % outputs (i.e. it represents a value of SIC % for which there is a 90 % probability that it is equal to or lower than the given SIC % level) versus OSI-430 SIC % observations for the 2015–2016 sea ice season. Some anomalies of various intensities can be detected in the series such as early 2015 melt-out events at Hall Beach and Churchill (▵), late 2016 freeze-up events at Akismi Island and Churchill (∘), and an early 2016 melt-out event at Hall Beach (◊).

The OFB site, near Iqaluit, has a behaviour that indicates an underlying error. In fact, neither the predictions nor the observations reach a SIC % value of 0 %, which is improbable in our study area. The reason behind these discrepancies with the validation data is that the OSI-409 product does not adequately estimate the SIC in this area, as can be seen in Fig. (7). The Frobisher Bay is usually ice-free around mid-July, which is never the case in the OSI-409 dataset. The source of this error is the land spill-over effect (i.e. land contamination) on estimated SIC combined with an inadequate sea ice presence estimation in the NSIDC sea ice monthly maximum extent climatology used as a mask for restricting areas where sea ice is likely to occur. Given the fact that the climatology mask states that there is possibly ice in the Frobisher Bay on that date, the algorithm attempts to measure it, with erroneous results. It has been found that this condition does occur in the Frobisher Bay and also west of Southampton Island, in Roes Welcome Sound.

.

Figure 7 Frobisher Bay OSI-409 error and time concordant MODIS True Color composite (Vermote, 2017)

Such errors make it important to emphasize that the results obtained from the model are to be used with care and ideally in combination with other sources of information such as local knowledge, other remote sensing imagery, and historical sea ice maps from national sea ice services.

4 Analysis of Hudson Bay sea ice spatio-temporal dynamics

The major asset of the IcePAC tool is that its output data give a probabilistic perspective on relevant sea ice event in comparison to the usual static descriptive statistics. Therefore, IcePAC gives not only the capacity to determine the mean event, but also to estimate the range of plausible events for a given site and date.

Here, the IcePAC outputs are used to assess the sea ice spatio-temporal dynamics, given different probability scenarios and in terms of three cover indicators, which are the length of the ice-free season (or its corollary, the ice-covered season), the probable complete melt-out week, and the probable complete freeze-up week.

4.1 Analysis with the IcePAC tool

Before presenting any results, the ice indicators analyzed in the next paragraphs must be clearly defined. First, the probable complete melt week corresponds to, for varying probability scenarios, the first week for which the SIC % is below 15 % in a given grid cell (starting its research from week 36 – 1 September – onwards). Second, the probable complete freeze-up week corresponds to, for varying probability scenarios, the first week for which the SIC % is above 15 % in a given grid cell (starting its research from week 14 – 1 April – onwards). To be considered valid, these events must be sustained for at least 3 consecutive weeks. Finally, the probable duration of the ice-free season corresponds to the gap, in weeks, between the different probable melt-out and freeze-up weeks.

The use of the 15 % limit to define the presence or absence of sea ice is a convention used by many authors (Andersen et al., 2006; Cavalieri et al., 1997, 1999; Divine and Dick, 2006; Gloersen et al., 1993; Pang et al., 2018) for SIC derived with passive microwave data and was therefore used in this analysis.

To estimate the aforementioned ice indicators, the probable SIC value for a given non-exceedance probability (p) was extracted from IcePAC for every location and week. For this analysis, a range of non-exceedance probabilities going from 5 % to 95 % was evaluated, with a step of 5 % between each analysis. Time series of the results were compiled and it is for these series that the complete freeze-up and melt-out events were identified. Figure 8 shows the estimated freeze-up and melt-out event weeks for p=50 %.

Figure 8 Estimated freeze-up (a) and melt-out (b) weeks for p=50 % in the HBS.

This detection process was repeated for each of the 20 738 grid cells of the simulation domain and for each probability step. It is worth noting that while the melt is described as a non-exceedance event (SIC < 15 %), directly deduced from p, the freeze-up is actually defined as an exceedance event (SIC > 15 %), deduced from 1−p.

Figure 9 Probable freeze-up and melt-out events for the coastal community of Puvirnituq (a), the range of probable ice-free season duration for the same community (b), and probable ice-free season lengths for all coastal communities of the HBS (c).

Figure 9a presents the probable freeze-up and melt-out events for the coastal community of Puvirnituq, located in the northeastern part of the Hudson Bay. In this figure we notice that melt has a 25 % probability to be completed by week 25 (18–24 June), only has a 10 % probability of being completed for week 24 (11–17 June), and is certain to be completed by week 31 (30 July–5 August). According to the two curves plotted in Fig. 9a, we can state that a complete freeze-up and complete melt at Puvirnituq can be expected, with a very high probability of occurrence (p>95 %) on weeks 50 (10–16 December) and 31 (30 July–5 August). Figure 9b shows an assessment of the range of probable ice-free season duration made for the same coastal community. Here, the freeze-up curve is inversed (1−P) as the shortest possible ice-free season duration is a combination of the latest possible melt (high exceedance probability) and the earliest possible freeze-up (high non-exceedance probability). By comparing the space between the two curves for the 5 % to 95 % probability range, we observe that the shortest possible ice-free season at Puvirnituq is 14 weeks and that the longest is 26 weeks. Figure 9c shows the ice-free season duration estimated for numerous coastal communities located in the study area using the method described for Fig. 9b. Particularly remarkable cases can be noticed, such as Cape Dorset, which displays a large variability in possible ice-free season duration given the shore lead polynya located in front of the community, Hall Beach in Foxe Basin with generally shorter ice-free seasons, with outliers corresponding to the occasional enlargement of the Hall Beach polynya, and finally Winisk, located in the south central part of the Hudson Bay, where ice tends to form early and melt late, which explains the shorter ice-free seasons when compared to other communities.

Table 1 Comparison between the CIS atlas and IcePAC p=50 % modelled occurrence weeks for the freeze-up and melt-out events at selected sites in the HBS. (“F” denotes freeze-up, “M” denotes melt-out and “NA” means not available.)

4.2 Comparison with the Canadian Ice Service atlas

The HBS has been analyzed since 1972 (i.e. first regional Hudson Bay map) by the Canadian Ice Service (CIS) and they have gathered numerous sea ice condition maps in the area from which they built a 30-year sea ice climatological atlas (CIS, 2013), portraying the average sea ice conditions and freeze-up and melt-out dates. As the information provided by the CIS atlas is given based on median SIC % values, a comparison with the p=50 % IcePAC output gives an outlook on the coherency of the model when compared with another source of data, built around a different methodology.

The “Dates of freeze-up and break-up” charts of the CIS 1981–2010 atlas depict the extent of ice on a biweekly basis during the freeze-up and break-up periods (CIS, 2013). They provide a pictorial representation of the evolution of the ice extent during those two periods. The freeze-up and break-up dates are estimated using median values for 1981–2010 computed from CIS regional ice charts, a collection of over 40 years of data which were digitized in the late 1990s as a raster with a 1 km grid size. These charts were prepared by trained sea ice analysts who used high-resolution RADARSAT images as inputs to their analyses, since its acquisitions started in 1996, in combination with AVHRR, NOAA, SSM/I, and ERS-1 data. However, all maps produced prior to 1996 were made without RADARSAT and therefore could be apprehended as less accurate than maps produced post-1996. The CIS also started using RADARSAT-2 in 2008. Another important information is that SIC % is estimated in tenths in CIS datasets in respect to the formalism of the egg code, meaning that SIC % marked as 1∕10 in an egg code could be in reality any value between 10 % and 19.99 %. Considering this, ice presence is defined as SIC $=>\mathrm{1}/\mathrm{10}$ for CIS data (CIS, 2013), while we stated that it was an SIC % $=>\mathrm{15}$ % for OSI-409.

Compared to the CIS data, IcePAC is based on algorithmically generated weekly averaged ice maps and calculated from adjusted frequency models. Also, there is a difference in spatial resolution between the two products; CIS uses 1 km digitized historical ice charts, while IcePAC uses passive-microwave-estimated SIC % resampled at a 12.5 km grid size. Finally, there is a time frame difference as the CIS atlas is built with data from 1981 to 2010 (30 years), while IcePAC is built with data from 1978 to 2015 (37 years).

The comparison, as displayed in Table 1, confirms that the freeze-up and melting dates identified by IcePAC are realistic when compared to the CIS historical data. Small differences in weeks are present and may be linked to a multitude of factors, the most important being the different methodologies used to generate the data.

The melt of the Hall Beach point does however come out as a relevant difference between the CIS and IcePAC melt week estimate. Since we find the freeze-up week adequately at the Hall Beach point and we use the same model distribution to derive both the freeze-up and melt information, it would be incorrect to simply flag this point as erroneous.

The overestimation of the ice-free season at Hall Beach by a 9-week gap compared to CIS could be explained by considering

• the land spill-over effect that would make the OSI-409 overestimate SIC %;

• the fact that the selected point (Hall Beach) is located on the edge of the polynya area;

• that the melt week statistics are not measured exactly for the same time period (1981–2010 versus 1978–2015);

• that we compare median values taken from CIS maps (values given in tenths of SIC %) with IcePAC outputs based on OSI-409 SIC % averages and that the definition of melt is not the same for both products (CIS is $<=\mathrm{1}/\mathrm{10}$, while IcePAC is <15 %).

By comparing the IcePAC P=50 % and P=15 % melt week output (Fig. 10) with the CIS melt dates, we can easily note that the CIS melt dates are quite variable in the polynya area and that by selecting a point a little further offshore, the results would have been comparable, towards advocation of the land spill-over effect as a coherent explanation. However, as we also compared the CIS melt weeks with a lower probability scenario (P=15 %), it turns out that the Hall Beach polynya does appear close to the community, as expected. This specific situation is certainly linked to the fact that IcePAC uses average SIC % values compared to the CIS median values. As the SIC % in this specific area tends to be either very low (open polynya) or very high (consolidated sea ice cover), it is credible to think that the average and median values do differ considerably. In a case of frequent low concentration like the Hall Beach polynya, the median tends toward a low SIC %, while the few higher SIC % events do bring the average SIC % up, increasing the gap between the median and the average values, enough so that with an average value we cannot detect the melt-out, therefore giving a plausible explanation for the 9-week difference noticed here.

Figure 10 Comparison of probable melt weeks for the coastal community of Hall Beach in Foxe Basin with the CIS melt weeks (a) and two IcePAC melt week scenarios for 50 % (b) and 15 % (c) probabilities. No data appear north of the IcePAC outputs since Hall Beach is located at the upper limit of the IcePAC model test domain.

4.3 Applicability of the IcePAC approach to other locations and data sources

The IcePAC modelling approach is replicable to any site for which SIC % data are available. Evidently, to ensure relevance for the resulting probability maps, the length of the time series used as inputs to adjust the beta models must be sufficient (∼30 years). Logically, the simulation domain grid cells will fall into one of the three following categories: (1) ice-covered grid cells with constant high SIC %, during the stable cover period; (2) marginal zone grid cells, in which SIC % will oscillate from high to low for the different years, during the freeze-up and melt periods; and (3) open water grid cells with constant low SIC %, during the ice-free season. As these are the ice regimes one could expect anywhere in the Arctic, including the HBS, there are no limitations on this side for using the IcePAC approach in other locations.

Other data sources could also be used with the IcePAC approach such as climate model outputs or different sea ice concentration maps. As climate model outputs provide future projections of SIC %, an evaluation of the range of probable SIC % patterns for a future year could be achieved (i.e. 2050 or 2080). However, one downside of these datasets is their coarse spatial resolutions.

5 Conclusion

The IcePAC tool permits an assessment of plausible sea-ice-related events for the entire range of probabilities for 20 738 sites (grid cells) in the Hudson Bay system, for all 52 weeks of the year. It is based on local (grid-cell) models that use the generalized beta distribution to describe the sea ice behaviour with four parameters at each site (with position and scale being fixed), based on historical 1978–2015 information from passive microwave imagery (OSI-409). From these parameters, IcePAC generates spatialized sea ice probabilistic information that can be used in any geographic information system or a web-based map interface for further analyses.

An analysis has been made to define, for each grid cell in the simulation domain, the plausible scenarios for each probability. A subsequent comparison with the 1981–2010 Canadian Ice Service sea ice climatology atlas (CIS, 2013) showed that the information generated with the IcePAC tool, for the p=50 % case, renders coherent probabilities for freeze-up and melt events over the HBS. A noticeable difference in the melt weeks was detected for the community of Hall Beach, suggesting that using median values instead of average SIC % could be of interest in order to be able to adequately detect specific events like polynyas. Another analysis, focused on the community of Puvirnituq, showed that it is possible to evaluate the range of plausible scenarios in terms of ice-free season length locally.

The model outputs generated with the IcePAC tool provide a novel probabilistic perspective regarding important events related to sea ice dynamics that was not available before. With its capacity to be utilized in other areas (with respect to the grid size of the passive microwave product), and the fact that it could easily be updated using new data, the IcePAC tool has the potential to provide valuable information on probable freeze-up and melt weeks as well as on ice presence and ice-free season lengths.

This relevant information will help decision makers such as engineers wanting to build a new marine coastal or offshore infrastructure to estimate ice hazards, fauna specialists trying to understand the vulnerability of a given species living dependent on the ice cover, or mariners wanting to estimate the feasibility of navigating a certain route. Turning the raw probabilistic information gathered from IcePAC into valuable thematic information will give stakeholders a capital gain in apprehending the possible spatio-temporal sea ice concentration patterns and in preparing for an effective mitigation of climate change impacts in coastal and offshore environments.