## 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.  
Figure 1: The NSR (red line) and representative sub-regions, as defined in the text, with sea ice conditions chosen from the ERA5 single-level monthly mean values for the date 1 July 2019.  
Figure 2: The grid point and statistical interpolations in sub-region No. 1.  
_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}{ll}Y_{1}&=\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}+\epsilon_{1}\\ Y_{2}&=\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}+\epsilon_{2}\\ \vdots\\ Y_{N}&=\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}+\epsilon_{N}\\ \end{array}\right., \tag{1}\]  
and in matrix notation:  
\[\left[\begin{array}{ll}Y_{1}\\ Y_{2}\\ \vdots\\ Y_{N}\end{array}\right]=\left[\begin{array}{