The results of computer simulation of the probability of accidents due to ship compressions by drifting ice along the Northern Sea Route (Eng)

Figure 1. Map of the route. Ice conditions are of the period 10-20 of May, 2018

In the model an accident occurs if vessel is in the zone of ice pressure and ice strength is upper than an amidship. If ridge concentration is more than 2 units than ice strength is taken to be equal to the greater value between level ice strength and strength of consolidated layer of ridges. During the model numerical simulations ice cover characteristics are set as stochastic variables (not deterministic), and their statistical distributions are model parameters.  The cumulative distribution functions are following: function of probability for vessel to get into zone under ice pressure; summarized length of a route in very close floating ice; relative lengths of parts of route in very close floating ice of various age; area of homogenous ice zones; lengths of parts of the route under ice pressure and without ice pressure; thicknesses of ice of various age; ridge concentration; ratio of ridge concentration to total ridge and hummock concentration; lengths of linear ridges; widths of ridges; ratio of ridge width to ridge altitude; thickness of consolidated layer of ridges. In Monte Carlo simulation, a random number generator provides a decimal fraction in the [0,1] interval which is considered as value of cumulative distribution function of the given characteristic of ice cover. Based on this value the quantile is estimated which is a specific value of ice characteristic. In the model we take into account only static interaction between hull and ice. Strength of the hull’s elements is estimated according to the requirements of the Russian Registry of Shipping. If randomly generated ice strength and its statistical distribution exceed the strength of a vessel, this cause hull destruction and emergency situation. Statistical distributions of ice characteristics are based on ice maps of the Arctic and Antarctic Research Institute (AARI) archive, data on AARI expeditions and literature sources. Ice compression strength depends on many factors.  Study on ice compression strength in ice interactions with different structures are high-demand and carried out by Russian [15-19] and foreign [20-24] specialists.

Ice strength depends on its thickness. During the period from the beginning of sustainable ice formation to the beginning of summer thawing ice thickness is determined by its age. Ridging develop predominantly in young and thin first-year ice (FYI). Newly-formed ridges are made of single ice rubbles that have the same thickness as the parent ice floe. Thus new ridges don’t pose a risk of ship besetting under ice pressure. The situation changes when the consolidation of part of single ice rubbles occurs and forms so called consolidated ice layer. It is supposed that navigators are able to avoid ridges if ridge concentration is up to 2 units according to Russian nomenclature (i.e. 2/5). In this case strength of level ice alone is taken into account; notably, both strength of deformation and strength of fracture are examined [25], and the lowest value is taken as the ice strength.

According to studies [26, 27] ridge concentration is 50% of total hummock and ridge concentration. Spatial orientation of ridges is random [27, 28], thus distribution of ridge orientations is set as proportional and mean angle between the direct course of a vessel and the ridge line is 45%.

Pre-processing of statistical distributions of summarized lengths of parts of the route characterized by specific parameters of ice cover is based on vector maps of ice conditions from AARI archive. These maps were made by the Centre of ice meteorological information of AARI on the base of remote sensing data [29].

Processing of vector map is made with use of ArcGIS as follows: polygonal objects of ice maps are crossed by linear object of navigation route. As a result a layer of linear objects is created that have the same sets of attributes as those of the crossed polygons. After that selection of objects according to their attribute values is done and lengths of the objects are calculated. Objects with total ice concentration not less than 9 units (9/10) with presence of ice older than young stage are identified, subsequently the selected objects are added to the new layer of parts of the route. Next step is to sort objects by age and to add them into separate layers. Lengths of individual objects and summarized lengths of objects that meet certain criteria are calculated. These result in ten-days series of lengths of the route segments in close ice, in close ice with presence of young (up to 10 cm thick), grey (10-15 cm thick), grey/white (15-30 cm thick), thin FYI (30-70 cm thick), medium FYI (70-120 cm thick), thick FYI (more than 120 cm thick). Besides, ten-days series of summarized lengths of parts of route in close ice with partial concentration of thick first-year ice of at least 5 tenths and summarized lengths of parts of route in close ice with sum of partial concentrations of thick first-year ice and medium first-year ice of at least 5 tenths are calculated.

The origin set of numeric variables for statistical distribution of parameter values of the model should be homogeneous. Thus the analysis of obtained numerical series for presence of interannual trend is made based on cumulative curves. Method of cumulative sums originally was applied in hydrology to examine the presence or absence of directed trends in the interannual dynamic of annual river discharges [30]. Afterwards it was applied for preliminary interannual variability analysis of any environmental parameters. Summary of this method is the following. A plot is made with years on the x-axis and cumulative sums of parameter values of specific years on the y-axis. In our case these are the values of the same within-year (intra-annual) ten-day interval of different years. Cumulative sum for particular year is the sum of values from the beginning of the time series to this particular year inclusive. For the first year it is the value itself, for the second year it is the sum of the values of both the first and the second year, etc. Data points of cumulative sums on the plot are connected with a line. The line must be close to the straight line in case of the absence of interannual variability. The kinked or broken-line curve show that there is a trend and the data is heterogeneous.

For curves that are close to straight line the original numerical series is divided into halves to examine its homogeneity. For kinked or broken-line curve the original numerical series is divided into parts in the points of inflection. If number of elements is not enough for statistical analysis, than original numerical series is divided into two parts in the point of the most significant inflection. Next step is testing the null hypothesis that two parts belong to one population, that means the absence of significant differences between these two samples. Verification is made in Mathcad with use of Mann-Whitney-Wilcoxon and Siegel-Tukey nonparametric tests.

Statistical distributions of model parameters were made in Mathcad based on methods of Hazen, Kritskiy-Menkel, Tchegodayev and Gringorten. Statistical distributions that were obtained by different methods varied inessentially. Therefore statistical distributions calculated by Gringorten method which combine three other methods were used in the model (Fig.2).

Figure 2. Cumulative distribution function of length of the route in close ice for the first ten-day period of May calculated by different methods.

Table 1 – Model outputs for different parameters. 

Thereby, parameters that governed the probability of accident are following: length of the route in close ice, total ridge and hummock concentration, probability of ship besetting under ice pressure and length of parts of route under ice pressure.

The results of 450 numerical simulations with varying vessel ice classes and vessel displacements are presented in Table 2.

Table 2 – Simulation results of navigations of different vessel ice classes and displacements

The results are visualized on Figure 3. The data indicates that stronger ice class reduces probability of accident. Thus, the model simulates correctly the increasing strength of the hull in case of increasing ice class.    

Figure 3. Simulation results for different ice classes of vessels 

Table 3 – Simulation results with forced setting of number of navigations 

 Conclusions. In this study the computer model of evaluation of accident probability due to ship besetting under ice pressure on the route “Sabetta ̶ Kara Strait ̶ Murmansk” was tested. To conduct numerical experiments the distributions of parameters of the model (i.e. characteristics of ice cover) are prepared, these are: thickness of consolidated layer of ridges, ratio of ridge width to ridge altitude, percentage of ridge concentration, total ridge and hummock concentration, lengths of parts of route without ice pressure, length of navigation route in close ice, relative length of route in close ice with presence of thick FYI, medium FYI, thin FYI and young ice. For data preparation ten-day interval vector maps of ice conditions on the route for the period 1997-2018 from AARI archive were processed. Obtained series of lengths were tested for having a trend based on method of cumulative curves and were tested for homogeneous by use of Mann-Whitney-Wilcoxon and Siegel-Tukey nonparametric tests. Statistical distributions were based on methods of Hazen, Kritskiy-Menkel, Tchegodayev and Gringorten.