Statistical Modeling of Arctic Sea Ice Concentrations for Northern Sea Route Shipping

Abstract

The safe and efficient navigation of ships traversing the Northern Sea Route demands accurate information regarding sea ice concentration. However, the sea ice concentration forecasts employed to support such navigation are often flawed. To address this challenge, this study advances a statistical interpolation method aimed at reducing errors arising from traditional interpolation approaches. Additionally, this study introduces an autoregressive integrated moving average model, derived from ERA5 reanalysis data, for short-term sea ice concentration forecasts along the Northern Sea Route. The validity of the model has been confirmed through comparison with ensemble experiments from the Coupling Model Intercomparison Project Phase 5, yielding reliable outcomes. The route availability is assessed on the basis of the sea ice concentration forecasts, indicating that the route will be available in the upcoming years. The proposed statistical models are also shown the capacity to facilitate effective management of Arctic shipping along the Northern Sea Route.

1. Introduction

Sea Ice Concentration (SIC) is a paramount environmental indicator of climatic change in the Arctic region. In light of the fluctuations in SIC, the Northern Sea Route (NSR) has become the focus of growing interest in recent years, as detailed in [1]. Consequently, accurate forecasts of SIC are crucial for the safe and efficient navigation of ships along the NSR, as demonstrated in studies by [2,3,4]. Chen et al. [5] evaluated the decadal changes in sea ice parameters and Arctic navigability, which suggested open water ships will be able to cross the Northern Sea Route and Northwest Passage between August and October during the period from 2045 to 2055, with a maximum navigable percentage in September. Liu et al. [6] calculated the navigational window of some Arctic routes based on sea ice density data from 2006 to 2015 and gained a better understanding of the sea ice conditions in the Northwest Passage, and the findings showed that the shortest navigational window for each route is 69 days. Zhou et al. [7] integrated satellite sea ice thickness data into the feasibility assessment model of the Arctic navigation, analyzing the navigability of the Arctic routes. The study showed that the overall navigational conditions of the Northwest Passage are poor, and the navigational window is relatively short. To reduce uncertainty in the prediction of the window period, Jun [8] researched the impact of radio detection observation on reducing uncertainty in Arctic sea ice prediction and identified issues regarding how to make forecast information convenient for crew and shipping companies to use. Li et al. [9] developed and verified two ship performance models which have been implemented in a voyage planning tool designed for summer Arctic operations of a cargo ship on the Northern Sea Route.

The SIC is typically recorded by satellite observations or by the observations of onboard master mariners. Researchers have utilized recorded data to formulate various models that attempt to discern the relationship between atmospheric anomalies and the sea ice regime. Stark et al. [10] incorporated the rheology mechanics of sea ice and its thickness distribution, ultimately positing that SIC is primarily governed by thermodynamics. Ye et al. [11] devised a scheme for reconstructing SIC using ice drift records, while Yu et al. [12] demonstrated that changes in SIC can be attributed to anomalous surface air temperatures and low-level wind fields. In a subsequent study [13], Ye et al. compared and appraised different methods for measuring sea ice type concentration. Mu et al. [14] proposed an ensemble-based Arctic Ice Ocean Prediction System, designed to facilitate forecasts of Arctic sea ice. However, physical models tend to be computationally intensive and require substantial amounts of data, and the resolution of available data may sometimes be insufficient to accurately model the geography along the NSR.

Furthermore, the statistical properties of sea ice concentration data are often neglected and under-investigated. As a statistical methodology, the autoregressive integrated moving average (ARIMA) model benefits from low computational requirements and high ease of use for time series analysis (TSA). Such techniques have been utilized in the context of oceanic environmental analysis, as demonstrated by Spanos [15], who compared an autoregressive (AR) model, a moving average (MA) model, and an autoregressive moving average (ARMA) model for simulating oceanic waves. Soares et al. [16] modeled significant wave heights through the use of an AR process that accounted for seasonality. Additionally, Agrawal and Deo [17] conducted a comparison between the ARMA model and an artificial neural network. These techniques have also been extended to a wide range of applications beyond oceanic waves, as seen in the works of Jiang et al. [18], who utilized the AR model for real-time prediction of ship motion, and Mao et al. [19], who applied the AR model for forecasting the speeds of container ships. The ARIMA model has been utilized in prior studies to examine the SIC. Piwowar and LeDrew [20] elucidated the autoregressive character of Arctic sea ice concentration. Ahn et al. [21] conducted a reanalysis of Arctic sea ice data through the ARIMA model, taking into account the changing temporal relationships between sea ice concentration and other climatic factors. Despite these studies, there remain limited instances where statistical models have been employed for the forecasting of Arctic sea ice concentration with a view toward the management of shipping along the Northern Sea Route.

This study proposes an ARIMA model to analyze and forecast the SIC along the NSR. A statistical interpolation approach is introduced to contract the SIC data obtained from the ERA5 [22] gridded reanalysis values into values along the NSR. The ARIMA model is trained and tested using reanalysis data from January 1980 to September 2022. Forecasts are generated for the period from 2022 to 2025 and validated through a comparison with ensemble experiments from the CMIP5 [23]. Based on this validation, an estimate for the 80% ice-free transit navigation window has been determined to serve as a guide for future NSR shipping. The results suggest that the proposed statistical methods are suitable and advantageous for the purpose of managing and planning NSR shipping.

This paper is organized as follows: Section 2 presents the statistical interpolation, and the ARIMA model is described in Section 3. In Section 4, the method is implemented as an example, and the results are discussed, while the conclusions are given in Section 5.

2. Interpolation of Sea Ice Concentrations

A typical NSR is depicted by the red trajectory in Figure 1. In accordance with the navigation plan of the ship that traversed this route, the NSR has been divided into eight segments, and consequently, eight sub-regions (depicted by green polygonal regions) have been established, each of which encompasses a single segment. In each sub-region, two longitudinal boundaries have been demarcated that define the starting and ending longitudes that the vessel traversed in accordance with the navigation plan; the northern boundary is situated 4 degrees north of the highest latitude reached on the route segment split in the section, while the southern boundary is located 4 degrees south of the lowest latitude reached on the route segment.

The transit navigation window, a critical determinant of a shipping company’s plans, is defined by the availability of the route. The acquisition of SIC data along the NSR from the ERA5 grid point values is indispensable in determining route availability. In the event that the grid resolution is adequate, interpolation from the nearest grid point values can yield acceptable SIC data along the route. However, such high resolution may not always be attainable, especially in areas of narrow straits. In such cases, interpolation necessitates a more comprehensive understanding of SIC in the entire sub-region, rather than being confined to the nearest grid point values. Additionally, when forecasting SICs for future shipping plans, the predicted high-resolution SICs often introduce significant uncertainties. In other words, it is easier to ascertain the presence of ice coverage along an extended route than to determine the SIC values at specific geographic locations. Consequently, the mean SIC of each route segment provides a more credible representation.

In Figure 2, the traditional method of interpolation via the nearest grid points is depicted on the left, while the proposed method, which leverages all grid values within the sub-region to enhance the interpolation, is presented on the right. A vast region provides significantly more information than a limited number of nearest grid points. For the purpose of statistical interpolation, the following features were selected from the gridded SIC values within a sub-region: total sea ice area, mean, variance, standard deviation, skewness, and kurtosis. A more extensive discourse on the statistical interpolation procedure is presented in the following section.

Statistical Interpolation of SIC

For each sub-region, let $Y_t$, $X_{t,sia}$, $X_{t,mean}$, $X_{t,var}$, $X_{t,std}$, $X_{t,skew}$, and $X_{t,kurtosis}$ denote the mean SIC of the route segment, area covered by sea ice, mean, variance, standard deviation, skewness, and kurtosis of SIC in the sub-region at time t, respectively. Thus, the mean SIC of the route segment is assumed to be described by the linear model as follows:

$$\left\{\begin{array}{cc}{Y}_1\hfill & ={\beta }_0+{\beta }_1{X}_{1,sia}+{\beta }_2{X}_{1,mean}+{\beta }_3{X}_{1,var}+{\beta }_4{X}_{1,std}+{\beta }_5{X}_{1,skew}+{\beta }_6{X}_{1,kurtosis}+{ϵ}_1\hfill \\ Y_2\hfill & ={\beta }_0+{\beta }_1{X}_{2,sia}+{\beta }_2{X}_{2,mean}+{\beta }_3{X}_{2,var}+{\beta }_4{X}_{2,std}+{\beta }_5{X}_{2,skew}+{\beta }_6{X}_{2,kurtosis}+{ϵ}_2\hfill \\ ⋮\hfill & \\ Y_N\hfill & ={\beta }_0+{\beta }_1{X}_{N,sia}+{\beta }_2{X}_{N,mean}+{\beta }_3{X}_{N,var}+{\beta }_4{X}_{N,std}+{\beta }_5{X}_{N,skew}+{\beta }_6{X}_{N,kurtosis}+{ϵ}_N\hfill \end{array}\right.,$$

(1)

and in matrix notation:

$$\left[\begin{array}{c}{Y}_1\\ Y_2\\ \phantom{⋮}\\ ⋮\\ \phantom{⋮}\\ Y_N\end{array}\right]=\left[\begin{array}{ccccccc}1& X_{1,sia}& X_{1,mean}& X_{1,var}& X_{1,std}& X_{1,skew}& X_{1,kurtosis}\\ 1& X_{2,sia}& X_{2,mean}& X_{2,var}& X_{2,std}& X_{2,skew}& X_{2,kurtosis}\\ \phantom{⋮}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \phantom{⋮}\\ 1& X_{N,sia}& X_{N,mean}& X_{N,var}& X_{N,std}& X_{N,skew}& X_{N,kurtosis}\end{array}\right]\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ {\beta }_2\\ {\beta }_3\\ {\beta }_4\\ {\beta }_5\\ {\beta }_6\end{array}\right]+\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ \phantom{⋮}\\ ⋮\\ \phantom{⋮}\\ {ϵ}_N\end{array}\right].$$

(2)

Finally, the linear model can be expressed in a compact form as

$$\mathit{Y}=\mathit{X}\mathbf{\beta }+\mathbf{\epsilon },$$

(3)

where $\mathit{\beta }$ represents coefficients to be estimated, $\mathit{X}$ is the design matrix, and $\mathit{\epsilon }$ denotes errors. $\mathit{\beta }$ can be estimated by minimizing the sum of squared residuals (RSS),

$$\begin{array}{cc}\hfill arg\underset{{\beta }_0,\cdots ,{\beta }_6}{min}RSS\left(\mathit{\beta }\right)& =arg\underset{{\beta }_0,\cdots ,{\beta }_6}{min}\left({\mathit{\epsilon }}^{\top }\mathit{\epsilon }\right)\hfill \\ \hfill & =arg\underset{{\beta }_0,\cdots ,{\beta }_6}{min}\left({(\mathit{Y}-\mathit{X}\mathit{\beta })}^{\top }(\mathit{Y}-\mathit{X}\mathit{\beta })\right)\hfill \\ \hfill & =arg\underset{{\beta }_0,\cdots ,{\beta }_6}{min}\left({\mathit{Y}}^{\top }\mathit{Y}-{\mathit{\beta }}^{\top }{\mathit{X}}^{\top }\mathit{Y}-{\mathit{Y}}^{\top }\mathit{X}\mathit{\beta }+{\mathit{\beta }}^{\top }{\mathit{X}}^{\top }\mathit{X}\mathit{\beta }\right),\hfill \end{array}$$

(4)

which leads to

$$\frac{\partial RSS\left(\widehat{\mathit{\beta }}\right)}{\partial \widehat{\mathit{\beta }}}=-2{\mathit{X}}^{\top }\mathit{Y}+2{\mathit{X}}^{\top }\mathit{X}\widehat{\mathit{\beta }}=0.$$

(5)

Finally, the estimator $\widehat{\mathit{\beta }}$ can be expressed as

$$\widehat{\mathit{\beta }}={\left({\mathit{X}}^{\top }\mathit{X}\right)}^{-1}{\mathit{X}}^{\top }\mathit{Y},$$

(6)

under the assumption that $rank\left(\mathit{X}\right)=k+1$ or

$$det\left({\mathit{X}}^{\top }\mathit{X}\right)\ne 0.$$

(7)

Additionally, strict exogeneity tells that the conditional expectation of $\mathit{\epsilon }$, given the design matrix $\mathit{X}$, is equal to zero, which implies that

$$\widehat{\mathit{Y}}=\mathbb{E}(\mathit{Y}\mid \mathit{X})=\mathit{X}\widehat{\mathit{\beta }}.$$

(8)

With Equation (8), a linear model has been developed to compute mean SICs of the NSR segments based on both the reanalysis SICs and their statistics.

3. Modeling of SIC and Time Series Analysis

The historical mean SICs of each NSR segment have been determined by the methods given in Section 2, forming a series of data points indexed in time order. This type of data can be described as a time series. In this section, we use the ARIMA model to analyze the SIC series.

3.1. Formulation of the ARIMA Model

The pre-requisite for a time series to be modeled by an ARIMA model is the demonstration of both stationarity and non-randomness. Stationarity greatly simplifies the complexity of the analysis and enhances the accuracy of estimation. Moreover, it helps to mitigate the risk of spurious regression [24]. However, a completely random series, such as white noise, is devoid of any correlations and therefore has no memory or influence on future values. That means that past information has no bearing on future information, and as such, it has no value for further analysis. Hence, before applying the ARIMA model, it is crucial to perform a stationarity test using the Augmented Dickey-Fuller (ADF) hypothesis test [25].

To represent the ARIMA model, we introduce the lag operator B first. When a SIC series y at time t is multiplied by B once, it means that time t has lagged one step. For instance,

$$\begin{array}{cc}\hfill y_{t-1}& =B\phantom{{}^1}{y}_t,\hfill \\ \hfill y_{t-2}& =B^2{y}_t,\hfill \\ \hfill & ⋮\hfill \\ \hfill y_{t-p}& =B^p{y}_t.\hfill \end{array}$$

(9)

A stationary and non-random time series contains the relevant information extracted by an ARMA model. The information at the current time t depends on the information in the past, which can be expressed mathematically as

$$\left\{\begin{array}{c}{y}_t=\underset{\mathrm{auto}-\mathrm{regressive}\left(\mathrm{AR}\right)}{\underbrace{{\phi }_1{y}_{t-1}+\cdots +{\phi }_p{y}_{t-p}}}+\underset{\mathrm{moving}\mathrm{average}\left(\mathrm{MA}\right)}{\underbrace{{ϵ}_t-{\vartheta }_1{ϵ}_{t-1}-\cdots -{\vartheta }_q{ϵ}_{t-q}}},\hfill \\ \phantom{y_t}={\phi }_{1}B{y}_t+\cdots +{\phi }_p{B}^p{y}_t+{ϵ}_t-{\vartheta }_{1}B{ϵ}_t-\cdots -{\vartheta }_q{B}^q{ϵ}_t,\hfill \\ {\phi }_p\ne 0,{\vartheta }_q\ne 0,\hfill \\ \mathbb{E}\left({ϵ}_t\right)=0,\mathbb{V}ar\left({ϵ}_t\right)={\sigma }_{ϵ}^2,\mathbb{E}\left({ϵ}_t{ϵ}_s\right)=0,s\ne t,\hfill \\ \mathbb{E}\left(y_s{ϵ}_t\right)=0,\forall s<t.\hfill \end{array}\right.$$

(10)

Equation (10) implies that

• the ARMA model can be written as ARMA $\left(p,q\right)$, which indicates that the highest order of the auto-regressive (AR) model is p, and the highest order of the moving average (MA) model is q;
• the error ${ϵ}_t$ is white noise; and
• the error at the current time is independent of past y.

If the series is differenced to obtain the stationarity, and the series contains seasonal components, the ARMA model can be extended to an auto-regressive integrated moving average model, which is expressed in a compact form as

$${\nabla }^d{\nabla }_S^D{y}_t=\frac{\Theta \left(B\right){\Theta }_S\left(B\right)}{\Phi \left(B\right){\Phi }_S\left(B\right)}{\epsilon }_t,$$

(11)

where

$$\begin{array}{cc}{\nabla }^d={(1-B)}^d,\hfill & {\nabla }_S^D={(1-B^S)}^D,\hfill \\ \Theta \left(B\right)=1-{\vartheta }_{1}B-\cdots -{\vartheta }_q{B}^q,\hfill & \Phi \left(B\right)=1-{\phi }_{1}B-\cdots -{\phi }_p{B}^p,\hfill \\ {\Theta }_S\left(B\right)=1-{\theta }_{1}B-\cdots -{\theta }_Q{B}^{QS},\hfill & {\Phi }_S\left(B\right)=1-{\varphi }_{1}B-\cdots -{\varphi }_P{B}^{PS}.\hfill \end{array}$$

(12)

The operator ${\nabla }^d$ is the ordinary difference component that represents the d-th differenced time series. The seasonal difference of order D is defined as ${\nabla }_S^D$, where S denotes the seasonal period, e.g., $S=12$ for a yearly seasonal period. $\Theta \left(B\right)$ and $\Phi \left(B\right)$ are the non-seasonal autoregressive and moving average lag polynomials, respectively. The remaining seasonal information is then extracted by the seasonal autoregressive and moving average lag polynomials, ${\Theta }_S\left(B\right)$ and ${\Phi }_S\left(B\right)$, respectively. Equation (11) is denoted as ARIMA $\left(p,d,q\right)\times \left(P,D,Q,S\right)$ in short.

3.2. Model Identification and Forecast

For an ARIMA $\left(p,d,q\right)\times \left(P,D,Q,S\right)$ model, the parameter d can be chosen by testing the stationarity of the d-th differenced series; S can be obtained from the seasonal decomposition, e.g., the Census X-11 method [26]. The remaining orders p, P, q, and Q can be determined from the properties of the AR$\left(p\right)$ and MA$\left(q\right)$ models, respectively.

Given an arbitrary AR$\left(p\right)$ model, from its partial autocorrelation function (PACF), we say that the AR$\left(p\right)$ model is “cutting off” at lag p. Consequently, p and P can be determined by this property from the PACF of series. Regarding an arbitrary MA$\left(q\right)$ model, its autocorrelation function (ACF) has the property of “cutting off” at q lags, which can help to find the order of q and Q. Once the orders have been determined, the parameters ${\vartheta }_i$, ${\theta }_i$, ${\phi }_i$, and ${\varphi }_i$ in Equation (11) are then estimated by using the maximum likelihood estimation (MLE) [27]. Finally, the minimum mean square error (MMSE) forecast [28] can be performed to produce a prediction.

4. Model Implementation and Validation

In this section, the statistical interpolation (described in Section 2) and the ARIMA model (described in Section 3) are applied and validated on data from ERA5 single-level monthly mean values for identifying the route availability. The workflow is illustrated in Figure 3. Sub-region No. 2 is chosen as an example and its context is described in detail.

4.1. Statistical Interpolation in Sub-Regions

The results of the feature importance analysis suggest that the SIC metrics of area, mean, variance and standard deviation are of significant relevance to the statistical interpolation, as depicted in Figure 4. Conversely, skewness and kurtosis do not seem to have any significant impact. Additionally, a close examination of Pearson’s coefficients, as illustrated in Figure 5, highlights strong linear correlations between the SIC’s area and mean, as well as between its standard deviation and variance. In this regard, the sea ice area, $X_{t,sia}$, and standard deviation, $X_{t,std}$, are selected out of the six statistics as features upon which to perform the interpolation. The results of the statistical regression are summarized in

Table 1. Statistical interpolation results.

	Coeff.	Std. Err.	t-Statistic	p-Value
${\widehat{\beta }}_1$	2.6803	0.056	47.970	0.000
${\widehat{\beta }}_2$	−0.2203	0.025	−8.899	0.000

, and the SIC value can be expressed as

$${\widehat{Y}}_t=2.6803X_{t,sia}-0.2203X_{t,std}.$$

(13)

In Table 1, the t-test suggests that the null hypothesis of $\mathit{\beta }=\mathit{0}$ can be rejected since p-values of ${\beta }_1$ and ${\beta }_2$ are equal to zero. Additionally, the $R^2$ value, which measures the proportion of variation in the dependent variables, is equal to $0.937$ and also verifies that the statistical regression fits the reanalysis data well.

Following the implementation of statistical interpolation across all sub-regions, the correlation between the interpolated and reanalysis SIC values is depicted in Figure 6. This illustrates that most of the points are situated in close proximity to the diagonal line, suggesting that the linear model effectively captures the SIC along the NSR. A diagnostic examination of the residual is depicted in Figure 7, which demonstrates, through the Q-Q plot and histogram, that the residual of the interpolation conforms to a normal distribution as intended. The scatter plots of “residual vs. fitted values” and “residual vs. $\mathit{X}$” evince randomness, thus affirming that the residual of the statistical interpolation constitutes white noise.

As previously noted in Section 2, the precision of the nearest grid point interpolation is strongly contingent upon the resolution of the grid. To evaluate the sensitivity of the statistical interpolation to grid size, the nearest grid point interpolation and statistical interpolation were performed with grid sizes ranging from 0.5${}^{\circ }$ by 0.5${}^{\circ }$ to 3.5${}^{\circ }$ by 3.5${}^{\circ }$. The nearest grid point interpolation with a resolution of 0.25${}^{\circ }$ by 0.25${}^{\circ }$ served as the reference. As depicted in Figure 8, it is evident that as the grid size increases, the nearest grid point interpolation progressively generates larger root-mean-squared errors (RMSEs, as indicated by Equation (14)).

$$\begin{array}{cc}\hfill & e_i=\sqrt{\frac{{\sum }_{t=1}^N{\left({\widehat{Y}}_{t,m,i}-{\widehat{Y}}_{t,0.25,1}\right)}^2}{N}},\hfill \\ \hfill & \mathrm{with}m\in \left\{0.5,1.0,1.5,2.0,2.5,3.0,3.5\right\},\hfill \\ \hfill & i=1,2,\hfill \end{array}$$

(14)

where ${\widehat{Y}}_t$ is the interpolated SIC defined in Equation (8). The nearest grid point interpolation and statistical interpolation are denoted as $i=1\mathrm{and}2$, respectively. The total number of timestamps is represented by N, and m represents the grid sizes. In fact, when the grid resolution changes from 0.5${}^{\circ }$ to 3.5${}^{\circ }$, the RMSEs increase by a factor of ten. However, the statistical interpolation exhibits a comparatively consistent RMSE trend, fluctuating around 0.05. While, at smaller grid sizes, the statistical interpolation yields higher errors than the nearest grid point interpolation. As the grid size exceeds 2 degrees, both interpolation methods produce comparable errors. When the grid size exceeds 3 degrees, the statistical interpolation proves to be significantly more accurate, yielding much lower errors in comparison to the nearest grid point interpolation. The insensitivity of the statistical interpolation to grid resolution, as demonstrated in Figure 8, expands its practicality and versatility, as it can be effectively employed with both high-resolution and low-resolution grids.

In this study, we undertake a comparison of the proposed statistical interpolation method with that of kriging interpolation with various variogram models. In order to establish a baseline for comparison, we obtain the average SIC values within each region through the use of the nearest grid point interpolation method, utilizing the highest resolution reanalysis of SIC data. Subsequently, the average SIC values are determined through statistical interpolation and kriging interpolation, and the root means square errors are calculated and compared against the reference values. The results (in

Table 2. The errors of statistical interpolation and kriging with different variogram models.

Region	Linear	Power	Gaussian	Spherical	Exponential	Hole-Effect	Statistical
1	0.139	0.110	0.101	0.117	0.093	0.076	0.043
2	0.178	0.143	0.182	0.199	0.199	0.157	0.060
3	0.071	0.076	0.100	0.101	0.081	0.102	0.036
4	0.221	0.202	0.206	0.196	0.201	0.134	0.035
5	0.157	0.154	0.159	0.162	0.154	0.171	0.022
6	0.178	0.181	0.189	0.159	0.151	0.135	0.040
7	0.124	0.118	0.199	0.152	0.130	0.143	0.029
8	0.307	0.222	0.262	0.276	0.239	0.135	0.048

) reveal that the average SICs obtained through the statistical interpolation method exhibit lower RMSEs in all regions as compared to the kriging interpolation.

Figure 9 depicts the October to December quarter mean reanalysis SIC (for reference), statistical interpolation, and kriging interpolation as an example. It is evident that the kriging interpolation produces the most pronounced deviation around the year 2010. Conversely, both the statistical and kriging interpolations produce results that are largely comparable, particularly after the year 2017. Hence, the proposed statistical interpolation can serve as an adequate alternative to the kriging interpolation. Furthermore, the statistical interpolation avoids the need for variogram model selection, yielding it more feasible to implement than the kriging interpolation.

4.2. Arima Modeling of SIC

The SICs along the given NSR obtained in Section 4.1 are plotted in Figure 10 as a time series. The time series contains a decreasing trend and a seasonal period of 12 months. To stabilize variances further, we preprocess SICs in two steps: (1) adding 1 to the original SICs and (2) performing a logarithmic transformation, i.e., $Y_t^{\prime }=Log\left(\widehat{Y_t}+1\right)$. The $Y_t^{\prime }$ is then used to establish the ARIMA model. Therefore, the original SICs can be obtained from the back-transformation, $\widehat{Y_t}=e^{{Y_t}^{\prime }}-1$.

Following the flowchart (Figure 3), the tests of randomness and stationarity are required next. From

Table 3. Ljung-Box white noise test of SICs in sub-region No. 2.

Lags	Ljung-Box Statistic	p-Value
1	286.935	0.0
6	807.854	0.0
12	1746.735	0.0

, the Ljung-Box test suggests that the preprocessed SIC series in sub-region No. 2 contains relationships at lags of 1, 6, and 12 months. In other words, the series is non-random.

Next, the Augmented Dickey-Fuller Test (ADF) tests the null hypothesis that a series is not stationary. Before carrying out the test, a seasonal difference, ${\nabla }_{12}^1{Y_t}^{\prime }$ (Equation (11)) is taken for the series. The results in

Table 4. Augmented Dickey-Fuller test of seasonally differenced SICs in sub-region No. 2.

Item	Value
p-value	0.019
test statistic	−3.208
critical value (1%)	−3.446
critical value (5%)	−2.868
critical value (10%)	−2.570

show that the p-value is smaller than 0.05, indicating that the seasonally differenced series are stationary. The statistic value is smaller than the critical values of 5% and 10%, and very close to 1%, which also supports the stationarity of the series.

In the ARIMA $\left(p,d,q\right)\times \left(P,D,Q,S\right)$ model, orders $d=0$, $D=1$, and $S=12$ have been confirmed. The orders p, q, P, and Q can be determined from the autocorrelation function (ACF) and the partial autocorrelation function (PACF) in Figure 11, the blue shaded regions represent the range of zero plus/minus two standard errors. When values fall within the shaded region, the ACF and PACF are considered to be zero. The lags of the red dots provide information to determine the orders of the ARIMA model. Within the first seasonal period, i.e., lag $h\in [0,12]$, the ACF cuts off at lags 3, 4, or 5, whereas the PACF cuts off at lag 2. In this regard, p is equal to 2, and q can be either 3, 4, or 5. Furthermore, the ACF has a peak at lag 12 and decays very rapidly at other multiples of the seasonal period, e.g., $h=12k$, for $k>1$. Therefore, Q is equal to 1. Nevertheless, the PACF decays slowly and oscillates at lags $h=12k$, for $k=1,2,\cdots$, and exhibits “tailing off” instead of “cutting off.” In this regard, P should be equal to 0.

Now only the value of q is absent in the ARIMA $\left(2,0,q\right)\times \left(0,1,1,12\right)$ model. To this end, the cross-validation is designed to find the best ARIMA model that gives the lowest Akaike information criterion (AIC, [29]). The total set of 492 months is split into a training set and a test set in five different ways. The size of the test size is always 60 months, while the training set sizes vary from 60 to 432 months (Figure 12), each of fold uses the same test set size but monotonically increasing training set sizes, the test sets cover the period from month 432 to month 492. The value $q=3$, which gives the lowest mean AIC among the cross-validations, is finally selected as the best candidate for the ARIMA model.

The model is then fitted with the in-sample data to compute the parameters in Equation (11). The size of the in-sample data set will influence the goodness of fit and the accuracy of the forecast. To choose a suitable in-sample size (i.e., training set size), we compare the RMSEs between the forecast SICs (obtained from training sets) and the true SICs (in test sets) in Figure 13, in each cross-validation, the model is first fitted by data in the training set. Then, the fitted model forecasts the SICs in the months covered by the test set. The RMSE represents the difference between the forecast and the true SICs in the same period. As the training set size increases, the RMSEs decrease monotonically. For training set sizes larger than 300 months, the RMSE curve shows convergence, following the recommended in-sample size in [30]. Eventually, 432 months are chosen as the in-sample size.

Using the MLE as described in Section 3.2 to solve for the parameters, the ARIMA $(2,0,3)\times (0,1,1,12)$ model is finally expressed as

$$\begin{array}{cc}\hfill (1-& B^{12})Y_t^{\prime }=\hfill \\ \hfill & \frac{1-0.5592B+0.5290B^2+0.3939B^3}{1-0.0851B-0.7502B^2}\left(1+0.8520B^{12}\right){ϵ}_t,\hfill \\ \hfill & {ϵ}_t\sim N\left(0,0.0028\right).\hfill \end{array}$$

(15)

Recalling the limitations of ${ϵ}_t$ in Equation (10), it is important to assess the residual of the ARIMA $(2,0,3)\times (0,1,1,12)$ model to ensure its validity. The residual is the difference between the actual SIC values and the values estimated by the model. The quarterly and yearly average residuals are shown in Figure 14, and they oscillate around zero, indicating a white noise pattern. The histogram, Q-Q plot, and autocorrelation function (ACF) of the residuals are plotted in Figure 15. The histogram and Q-Q plot demonstrates that the residuals are almost normal in distribution, and the ACF indicates that the residuals are independent of time, which are the two main assumptions for the residual in the ARIMA model. Hence, the ARIMA $(2,0,3)\times (0,1,1,12)$ model has effectively captured most of the information from the original SIC series and the residual is, in essence, white noise.

Since the ARIMA model has demonstrated acceptable goodness-of-fit, we employ the MMSE forecast [28] to project the SIC from 2021 to 2025 (in Figure 16). The red lines depict the input data obtained from reanalysis, while the blue dotted lines serve as a reference for verification. The green lines represent the SIC values estimated via the ARIMA model, and the purple lines indicate the forecasts made by the same model. The red shaded areas denote the 95% confidence intervals of the forecasts. Prior to 2021, the estimated SIC values are nearly identical to the input data from reanalysis. Additionally, the forecasted SIC values for the year 2021 align well with the reanalysis data, signifying that the proposed ARIMA model is capable of generating reliable predictions of SIC.

Furthermore, the quarterly means of the estimated and forecasted SIC values are illustrated in Figure 17. The red lines denote the reanalysis data, while the blue curves represent the SIC values from the ARIMA model. The red shaded areas represent the 95% confidence intervals of forecasts. The plotting of the quarterly means of SIC values allows for easier discernment of the changing trends of the SIC across the quarters. Throughout all quarters, the estimated SIC values align with the reanalysis data, effectively capturing the declining trend of SIC values. The largest discrepancy occurs in the October to December quarter (OND) of the year 2021. Nonetheless, the reanalysis SIC remains within the confines of the 95% confidence intervals.

The distributions of quarterly means of SIC values across the various sub-regions are depicted in Figure 18. It is observed that the SICs in sub-region No. 1 are notably lower than in other sub-regions, reaching close to zero from October to December. Conversely, higher SICs are observed in the other sub-regions from January to March and April to June. Additionally, a discernible downward trend in SICs is evident from July to September. The forecasting of the SIC values reveals that the fluctuations are less pronounced than in the reanalysis data, yet the trends are properly captured.

4.3. Validation of SIC Forecasts

We have shown the high goodness-of-fit of the proposed ARIMA model. Nevertheless, the forecast SICs need further validation. To this end, SICs are selected from nine CMIP5 experiments [23]. Although the CMIP5 lacks data before 2006, both CMIP5 and ERA5 provide SICs from January 2006 to September 2022. Therefore, this period is chosen to compare the accuracy of the proposed model. Figure 19 identifies the absolute errors of the monthly averages of SICs between CMIP5 and ERA5. Each dot represents an absolute error in a single month, and the color represents the quarter to which the month belongs. Figure 20 compares the RMSE of the monthly average of SICs between CMIP5 and ERA5. The most accurate experiment is listed at the top and the least accurate experiment at the bottom. The top four experiments are selected as the reference to validate the forecast SICs of the proposed ARIMA model.

In October, November, and December in sub-region No. 2, the monthly averaged SICs of the selected experiments are plotted in Figure 21. The solid blue line represents the mean SICs of these four experiments, surrounded by the 95% confidence interval (green shaded region). The ERA5 SICs are plotted as the solid brown line, and the ARIMA model regressed SICs are represented as the dashed blue line. From 2006 to 2015, the mean experimental SIC deviates from the SICs of ERA5 and the ARIMA model. Between 2015 and 2022, all three lines lie close together. From 2023, the ARIMA model forecast SICs stay close to the mean SICs of the CMIP5 experiments, the differences of which are bounded by 0.1 of the SIC value. Therefore, we can conclude that the proposed ARIMA model can provide the same confidence for forecasts as CMIP5 experiments.

4.4. Estimation of the Transit Navigation Window

The reanalysis and ARIMA-model-forecasted SICs serve as the basis for determining the navigational window for transit. SIC values that fall below 0.15 or 0.25 are considered to mark the threshold of ice-free conditions, as suggested in [31]. Hence, an 80% navigational window for transit is defined as the number of days per year in which the SIC values along 80% of the route surpass the threshold of ice-free. For any given date, if the SICs in six of the eight sub-regions fall below the threshold, that date is deemed to be included within the navigational window.

Pursuant to this principle, we proceed to estimate the 80% transit navigation window from the year 2022 to 2025. Figure 22 portrays the reanalysis and forecast transit navigation windows in conjunction. The navigation window that results from the use of a higher ice-free threshold is wider than the estimate obtained with a lower threshold, though the trends of both navigation windows are consistent. The navigation windows undergo cyclical fluctuations from wide to narrow from the year 1991 to 2025. Prior to 2013, the navigation window would experience contraction approximately every seven years, in 1999, 2007, and 2013. Since then, the navigation window has remained periodically vibrational, yet the narrowest windows are now noticeably broader than in the past. In 2021, the navigation window commenced in early August and terminated in mid-September (when using a high SIC threshold), implying a relatively narrow window. By the year 2022, the commencement of the window was postponed by approximately two weeks, but the window was extended until October. Since 2023, the navigation window will broaden, remaining from mid-July to the end of September. The navigation window will continue to expand slightly until 2025.

To validate the estimated navigation window using SICs from the statistical interpolation, we also utilized the SICs from the kriging interpolation to assess the navigation windows. As demonstrated in Figure 23, the navigation windows obtained from both kriging and statistical interpolations are in agreement for the period from 2016 to 2021. This convergence attests to the suitability and veracity of the proposed methodology for estimating the transit navigation window.

5. Conclusions

This study introduces a systematic method for forecasting Sea Ice Concentrations along the Northern Sea Route for Arctic shipping and evaluating the transit navigation window. The Arctic region is divided into sub-regions that encompass the NSR segments, and the mean SIC of each NSR segment in each sub-region is obtained through the proposed statistical interpolation. This technique is found to be less sensitive to grid size variations than traditional nearest grid point interpolation, and yields stable results even with low grid resolution, as it incorporates all statistical information within the sub-region. Moreover, this proposed statistical interpolation has been demonstrated to exhibit comparable accuracy to the widely employed kriging interpolation method.

The interpolated SICs are subsequently analyzed and modeled through the use of an Auto Regressive Integrated Moving Average model, which proves capable of producing accurate short-term SIC forecasts. This purely statistical model, which only relies on reanalysis (historical) data, has a forecast accuracy comparable to that of the more complex Climate Model Intercomparison Project 5 physical models, which require substantial meteorological and oceanographic data. Moreover, the ARIMA model is relatively simple to understand and implement, and is computationally efficient, making it an ideal candidate for real-time SIC estimation systems in onshore and offshore environments. The ARIMA model can be updated in real-time by re-training it with new SIC data, allowing for dynamic correction of forecasts, and continuous creation of new forecasts as the training set expands.

Finally, the SIC forecasts are applied to evaluate the transit navigation window for the years to come, indicating that the selected NSR is poised to become readily accessible to shipping, with improved sea ice conditions in the years ahead. In this investigation, we selected a grid size of 0.25${}^{\circ }$ by 0.25${}^{\circ }$ in order to assess the transit navigation window and authenticate the proposed statistical models. Although this grid resolution is relatively large and may augment the uncertainty in estimating the transit navigation window, the utilization of high-resolution sea ice information through the proposed statistical technique can yield a more precise transit navigation window. The comparison of uncertainty between high and low resolution represents an intriguing topic and leaves further examination.

The present study proposes a methodology that presents a reliable and efficient approach to forecasting sea ice conditions along the Northern Sea Route. The suggested technique is characterized by its explicitness and low computational requirements, which enable its successful implementation both onboard and onshore. The auto-regression model’s inherent characteristics facilitate continual learning and updating, thereby providing a sound basis for weather forecast justifications. The proposed methodology has significant implications for shipping companies in terms of managing and scheduling their Arctic shipping itineraries. By providing a preliminary estimate of the navigation window, the methodology allows companies to optimize their shipping operations along the NSR, resulting in reduced risks of delay and financial loss.