CC BY-NC-ND
0ELSEVIER
Contents lists available at ScienceDirect
Cold Regions Science and Technology
journal homepage: www.elsevier.com/locate/coldregions
Estimating operability of ships in ridged ice fields
Lauri Kuulialaa'*, Pentti Kujalaa, Mikko Suominena, Jakub Montewkaa-b-c
aAalto University, School of Engineering, Department of Mechanical Engineering, P.O. box 14100, Aalto FI-00076, Finland b Finnish Geospatial Research Institute, Department of Navigation and Positioning, Masala 02431, Finland cGdynia Maritime University, Faculty of Navigation, Department of Transport and Logistics, Gdynia 81-225, Poland
ARTICLE INFO
ABSTRACT
Article history:
Received 27 March 2016
Received in revised form 17 November 2016
Accepted 7 December 2016
Available online 12 December 2016
Keywords: Winter navigation Ridged ice Ship's performance Transit simulation Statistical simulation Sensitivity analysis
A method for estimating ship's resistance caused by sea ice ridge keels is revised and used as a part of a method for predicting performance of ships in ridged ice conditions. The resistance method is based on a continuum plasticity model of ridge rubble and is simple to compute. The performance prediction method combines deterministic simulations of ship motion with probabilistic modelling of ridged ice fields. Performance estimates given by the model are distribution of attainable mean speeds for given ice conditions and probability of the ship being able to operate independently.
A comprehensive sensitivity analysis was performed to gain insight into the model and identify possible problematic parameters. The sensitivity analysis covered both the ice conditions and modelling assumptions.
Two data-sets were used to test the simulation method. One set included the depth profile of sea ice, machinery data and the speed of a ship operating in ridged ice. The resistance method was able to predict the mean speed over 3km well. The second data-set consisted of a history of ship's speeds and positions from AIS data and ice conditions estimated by a numerical ice model HELMI, developed in the Finnish Meteorological Institute. Observed mean speeds were mostly well within the distributions of mean speeds simulated by the transit simulation model. Predictions of independent operation were also promising.
© 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND
license (http://creativecommons.org/licenses/by/4.0/).
1. Introduction
Ridged ice is a prevalent form of sea ice. In the Baltic, up to one third of the total ice volume is contained in ridges with maximum depths varying typically between 5-15 m. (Lepparanta et al., 1995) Ice ridges form in areas of local compression of the ice field, typically when leads or channels are closed by compressive forces caused by wind and currents.
Sea ice ridges constitute a major hindrance to vessel traffic in ice-covered waters and ridged ice is among the most difficult conditions encountered by ships navigating in ice. In the Baltic, there is a winter navigation system in operation into Finnish and Swedish ports. Meanwhile, Arctic waters are experiencing increased activity due to the shrinking of multi-year sea ice cover. The Northern Sea Route between Europe and East-Asia along the northern coast of Russia is expected to become an increasingly important transport route. The ability to predict ship performance in ridged ice is important for the planning of operations and transport routes.
* Corresponding author. E-mail address: lauri.kuuliala@aalto.fl (L. Kuuliala).
Ridging starts as rafting and continues as such until the stresses in the ice sheets, caused by friction and roughness of the ice cover, fail the ice in buckling and subsequently in bending. Pile-ups of broken ice pieces are formed under and over the parent ice sheet. Ridging continues until the driving forces can no longer sustain it or the strength of the parent ice sheet is exceeded elsewhere, resulting in a new ridging event. (Tuhkuri et al., 1999)
The underwater part of the ridge is called the keel and it consists of broken ice pieces and possibly a layer of rafted ice near the waterline. In mature ridges there is a consolidated layer in the keel, which extends down from the waterline. In the consolidated layer the water filling the voids between ice-blocks is frozen. The consolidated layer is inhomogeneous and it has a varying thickness that is on average 1.5-1.7 times the surrounding level ice thickness. (Bailey et al., 2015; Tuhkuri et al., 1999)
The ice pile-up on the surface is called the sail and it is composed of ice blocks that may be frozen together at the contact areas. The voids between blocks are filled with air or snow. (Strub-Klein and Sudom, 2012)
Sail and keel masses are in hydrostatic balance over large length scales but locally significant deviations from isostatic balance have
http://dx.doi.org/10.1016/jxoldregions.2016.12.003
0165-232X/ © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by/4.0/).
been observed (Kankaanpaa, 1989; Tin and Jeffries, 2003; Veitch et al., 1991). The ice cover surrounding ridges can experience deformations and carry significant stresses due to local hydrostatic imbalances (Bowen and Topham, 1995). Due to the role of rafting in ridge formation, there can be large lateral displacements between sails and keels (Melling et al., 1993).
Ship's resistance in ice ridges arises mainly from breaking the consolidated layer and displacing the keel rubble. The performance of a ship in ridged ice fields is influenced by the sizes and spac-ings of ridge keels. The spacings of ridge keels can be modelled by the exponential (Hibler et al., 1972) or the log-normal distribution (Wadhams and Davy, 1986). These models assume ridge spacings to be independently distributed. Lensu (2003) proposed a ridge spacing model, where ridge spacings are correlated and clustering of ridges is taken into account. Ridge sail heights can be modelled by the exponential distribution (Lepparanta, 1981; Lewis et al., 1993; Wadhams and Horne, 1980). The keel sizes and spacings are not directly observable in operational situations, causing significant uncertainty in predicting the performance of a ship in ridged ice fields.
Methods for predicting ship's resistance in ice ridges have been presented by Keinonen (1979), Mellor (1980), Malmberg (1983), Abdelnour et al. (1991) and Riska et al. (1997). The method by Malmberg (1983) has been used in the transit simulation models (La Prairie et al., 1995; Patey and Riska, 1999) and the method by Riska et al. (1997) in route optimization models (Guinnes et al., 2014; Kotovirta et al., 2009). Kamesaki et al. (1999) used the simulation model given by Patey and Riska (1999) in a transit simulation of the Northern Sea Route, that covered also economic aspects of the studied routes. Mulerhin et al. (1999) presented a transit simulation model, that takes ridged ice into account. They used look-up tables made using expert input and various data sources to assess ships' operation in ridged ice. Montewka et al. (2015) developed a data-driven probabilistic model for predicting ship performance in ridged ice by Bayesian networks.
Because of the amount of uncertainty concerning the ridge fields that a ship encounters, it is reasonable to model the performance of a ship as a stochastic phenomenon. The transit simulation model by Riska et al. (1997) is completely deterministic and while the method by La Prairie et al. (1995) generates the ridge fields for the simulation using statistical methods, the result does not account for the variation of results arising from the differences in realized ridge field geometries. Patey and Riska (1999) generated a large number of long transits for each investigated ridging condition to make an estimate of the average speed in the given conditions. The analysis was explicitly limited to ice conditions where the ship was likely to be able to operate independently, thus limiting severe ridging outside of the analysis. Kamesaki et al. (1999) extended the analysis to model distributions of mean speeds in ridged ice.
In this article, a similar method to those used by La Prairie et al. (1995) and Patey and Riska (1999) is extended to produce a statistical description of ship performance in ridged ice conditions. Performance is defined here as attainable mean speed and possibility of unassisted operation in given ice conditions. Mean speed is not given as a point estimate as in (La Prairie et al., 1995; Patey and Riska, 1999) but as a distribution reflecting the uncertainty of describing ridged ice with a small number of parameters similarly to (Kamesaki et al., 1999). Further, the analysis is not limited to conditions, where the ship is likely to operate independently but rather an estimate of the probability of becoming beset in ice is provided. In the modelling approach adopted in this paper, the length of the simulated leg and initial speed of the ship become also significant factors. These results are considered to be novel. Finally, a mapping from ice forecast data to ridge spacing and size distributions is proposed.
A comprehensive sensitivity analysis of the model is also presented. This expands the work presented in (Patey and Riska, 1999)
by extending the sensitivity analysis from the ice conditions also to the assumptions of the model. As far as the authors are aware, a sensitivity analysis of this extent has not been published for a similar model for ship's performance in ridged ice.
Section 2 describes the modelling approach used for the resistance of the ship, ridged ice and performance of the ship. An overview of the data-sets used in testing the model is given in Section 3 and the sensitivity analysis is outlined in Section 4. The results are presented in Section 5. The modelling assumptions and results are discussed in Section 6 while Section 7 contains a summary of the work.
2. Modelling ship transit in a ridged ice field
2.1. The equations of motion
Operating of a ship in ice is modelled by solving one degree of freedom equations of motion of a rigid body. Only the surge direction is considered. Excitation is composed of resistance and net thrust. Resistance is divided into components caused by breaking of level ice and the consolidated layer, displacing ridge keel mass at the bow and the sliding of keel mass along the parallel midship. The effect of ridge sails to the resistance of the ship is disregarded as it is assumed to be small compared to the overall resistance. It is assumed that the power setting, propeller revolutions and pitch remain constant when the ship is transiting ridged ice. The ship is further assumed to follow a straight line without trying to optimize its course over ridge fields. These significant simplifications are introduced to simplify modelling.
The ship is modelled with a set of parameters describing its main dimensions, bow shape and thrust. The parameters and their symbols are given in Table 1 and illustrated in Fig. 1. For the simulations, the lengths and angles were measured from lines drawings of the ships. Waterline opening angle, a, is the angle between the tangent of the waterline and mid-line of the ship at the location of the bow, where the breadth is half of the maximum breadth. The stem angle, 0, is measured at the waterline from the profile view. Normal angle at the horizontal at the bow, X, is calculated by X = arctan(tan0/sina). Coordinate x tracks the position of the ship's bow along the route of the ship.
Table 1
Parameters used to model the ship.
Parameter Symbol Unit
Length of waterline L m
Length of bow entrance Lb m
Length of parallel midship Lm m
Breadth at waterline B m
Draught T m
Waterline half-angle a rad
Stem angle 0 rad
Normal angle to horizontal at bow x rad
Mass of displacement m kg
Bollard pull Tb N
Open water speed vow m/s
Numerical integration is performed by Newmark's method. Position and velocity for time step j are given by
xj = Xj_1 + Vj_1dt + Q aj_1 + 1 j dt2
Vj = V-1 + 2(aj-1 + aj)dt
where dt is the length of time step and Vj and aj are velocity and acceleration at time step j. (Newmark, 1959) Acceleration at time step j is
aj = m,
where Fj is the excitation at time step j and m is the mass of the ship. Hydrodynamic added mass and damping are not taken into account in this model as they are assumed to be small compared to the other excitation components. Because the resistance components are dependent on velocity and position, the acceleration, velocity and position of a time step are interdependent. Iteration has to be carried out within each time step to solve the equations of motion. The first round of iteration is calculated assuming that a° = aj-1, where the superscript indicates the iteration round. In further rounds of iteration, the excitation is solved at the velocity and position calculated with the acceleration of the previous iteration round. Iteration is stopped when the relative change in the acceleration between rounds of iteration is less than 0.001.
Excitation at time step j is calculated by
Fj(Xj, Vj) = Tn(Vj) + Rb(Xj) + Rm(Xj) + Rl(x, Vj),
where Tn is net thrust, Rb and Rm are the ridge resistance components caused at the bow and midship respectively and Rl is the resistance caused by breaking level ice.
The ridge resistance components are calculated by the method presented by Malmberg (1983). The method is derived by using Rankine's plasticity model for the ridge rubble, which is assumed to be a cohesionless continuum material. The internal friction angle, 0r, is needed to characterize the plasticity behaviour of the material.
The resistance component due to displacing the ridge keel at the bow is calculated by
Rb =Cphr(x)(0.5B + hr(x)tanX cos a) x (u cos a + sin X sin a),
The integration is performed over the parallel midship area from the bow shoulder to the stern shoulder and the term (hr(x)/T - 0.5)B is analysed only when hr(x) > 0.5T. Cm is a constant dependent on the Poisson's ratio, m, of the ridge rubble given by
Cm = Ui(1 - P)PDg
1 - V'
is the coefficient of friction between ship and ice.
According to Mellor (1980), m can have values between 0.21-0.3. Cm has values in the range of 35-77 kg/m2 s2. The value of Cm is set to 45.9 kg/m2 s2 for the simulations (La Prairie et al., 1995).
The resistance components due to breaking level ice is calculated using the method by Riska et al. (1997). The level ice resistance is assumed to be linearly dependent on the speed of the ship:
Rl = C1 + C2Vj,
where C1 and C2 are given by
C1 =/1 (2T + 1 BLmhj(x) + (1+ 0.021*)
x (/2Bhi(x)2 + /3Lbhi(x)2 + /4BLbhj(x))
C2 =(1 + 0.063*)(g1hj(x)15 + g2Bhj(x)) 1.2T) B2
+ g3 1 +
where hi(x) is the ice thickness at the bow. The same resistance formula is also used for the resistance caused by the consolidated layer, which is assumed to have a thickness of 1.5 times the level ice thickness. Coefficients f1 .. .f4 and g1 ... g3 are given in Table 2.
The thrust of the ship is modelled using net thrust, i.e. the difference of thrust and open water resistance at a given speed. Net thrust is a function of ship speed and it is calculated by
j *( < - - 2 ( VW )}
where Vj is the speed at time step j, vow is the open water speed of the ship and Tb is the bollard pull, which can be estimated by
(5) Tb = Ke(PdDp)2'3,
where hr(x) is the local ridge thickness at the bow. Cp is a constant that depends on the internal friction angle of the ridge rubble calculated by
Cp = (1 - p)p&g
1 + sin 1 - sin '
pA is the difference of the densities of sea water and ice.
Porosity of ridges varies between 0.25-0.4 and the internal friction angle between 47°-58°. Density of sea-water ice is in the range 890-930 kg/m3 and the density of sea-water is about 1025 kg/m3. (Mellor, 1980) The value of Cp varies within 4000-11000 kg/m2 s2 with this range of parameter values. Following La Prairie et al. (1995), the value of Cp is set to 7500 kg/m2 s2 for the simulations.
The midship resistance component is given by
Rm = CmT I
i>x)+( T - 0-5) B)
where Pd is power in kW, Dp the propeller diameter in m and bollard pull is given in kN. Ke is a coefficient that takes the number and type of propellers into account. Values of Ke for merchant vessels are given in Table 3. (Riska et al., 1997) Net thrust varies from zero at open water speed to bollard pull at zero speed. The interaction of the propeller and ridge rubble decreases propulsive efficiency, which is modelled by multiplying net thrust with a coefficient Ctf (Kostilainen, 1981; Leiviska, 2004). The coefficient can have values between 0 and
Table 2
Coefficients of the level ice resistance formulas (Riska et al., 1997).
/1 230 N/m3
/2 4580 N/m3
/3 1470 N/m3
/4 290 N/m3
g1 1890 N/(m2-5/s)
g2 670 N/(m3/s)
g3 1550 N/(m3-5/s)
Table 3
Values of Ke for estimating bollard pull. (Riska et al., 1997).
Number of propellers Fixed pitch Controllable pitch
1 0.7 0.78
2 0.88 0.98
3 1.01 1.12
1 and it is initially set to 1. Other values of the coefficient are included in the sensitivity analysis.
In conditions where the ship can not maintain forward speed, ramming of ridges is modelled using the general approach described by La Prairie et al. (1995). When the speed of the ship drops to zero, it reverses thrust and backs up for a pre-determined distance. The distance is currently set to two ship lengths. During reversing, the level ice component is not calculated as it is assumed that the ship is reversing in the broken channel. The power setting is full astern while the ship reverses (La Prairie et al., 1995). Reversing causes a decrease in propulsive efficiency, which is modelled by multiplying the net thrust astern by a coefficient Ctb, which is the ratio of thrust astern to thrust forward at full power. Ctb is initially set to 1, but lower values are included in the sensitivity analysis.
After the ship has reversed for the determined distance, it sets the thrust forward with full power and starts to accelerate in the pre-broken channel towards the ridge. The level ice resistance is not calculated until the ship reaches the end of the broken channel. If the ship does not have sufficient thrust to reverse after coming to a stop, it is considered to be beset in ice.
The sensitivity of the method to the length of time step was tested. The changes in the predicted mean speed over 15 minutes was negligible for dt < 0.15 s. The time step was set to 0.1 s for the simulations.
2.2. Modelling of ridged icefields
In this section, the method used for modelling ridged ice fields is presented and the simplifications of the model are identified. The parametrization of the model using obtainable forecast data is also outlined. Finally, the method for making estimates of ship performance in ridged ice fields is presented. The model is limited only to full ice coverage as performance of ships in ridged ice conditions is studied.
In the most simple form, ridges are modelled as having triangular sails and keels. The model characterizes a ridge with sail height hs, keel depth hk, slope angles for sail and keel, Ks and Kk respectively as well as porosity, p. The ratio of keel depth to sail height for first-year sea ice ridges is 4-5 based on hydrostatic balance and assumptions of keel and sail angles (Wright et al., 1978). A large study of data on first-year sea ice ridges gave values for the hk/hs ratio of 4.35-5.2 (Strub-Klein and Sudom, 2012). Parameters of the ridge model are shown in Fig. 2.
Porosity of ridge keels varies between 0.25-0.4. Typical values for the keel slope angle are reported to be 25-30°. (Mellor, 1980; Timco and Burden, 1997). Strub-Klein and Sudom (2012) give the average ratio of keel depth to keel width to be 4.85, which leads to Kk = 22° assuming a triangular cross section of the keel. The cross-sectional area of a ridge keel is
Ak = V hk tan Kk.
Depths of ridge keels are limited by the increasing energy needed to submerge the ridge keel further. When the parent ice sheet is not strong enough to penetrate into the forming ridge, it will break in front of the ridge resulting in lateral growth of the keel or the start of another ridge further away. The limit depth is mainly determined by
Fig. 2. Structural model of ridges used in generating random ridge fields for the model. hi and hc are the thickness of level ice and consolidated layer. js and jk are the sail and keel opening angles while hs and hk are the height of sail and depth of keel.
the thickness of the parent ice sheet (Hopkins et al., 1991; Parmerter and Coon, 1972). The limiting depth of ridge keels for the Baltic is given by Kankaanpaa (1997) to be
hlim = 17.64h°5,
where hp is the thickness of the parent ice sheet.
Sail heights and keel depths of sea ice ridges can be modelled to be distributed according to the exponential distribution. The distribution was first proposed for the Arctic (Wadhams, 1980) and has been confirmed also in the Baltic for sail heights (Lepparanta, 1981; Lewis et al., 1993). The height distribution is parametrized by a cutoff height, hc and the shape factor of the distribution, k. The mean height of sails higher than hc is hc + kand standard deviation of heights is k-1. The cut-off height is needed to distinguish ridges from the noise of the measurement system used to observe sail heights and keel depths for determination of k. The cut-off height for sails in the Baltic is typically 0.4 m and kis 0.1-0.2 (Lepparanta, 2005). The distribution for ridge sail heights greater than the cut-off size is
p(h; hc, k) = k exp(-k(h - hc)).
It is assumed, that ridge keel depth is on average related to sail height and that the sail height distribution can be used to model keel depths. When ridge fields are generated for simulations, the sail heights are multiplied by the ratio of keel depth to sail height to obtain the keel depths. Value of 5 is used for the keel depth to sail height ratio.
Ridge keel spacings can be modelled by the exponential distribution (Hibler et al., 1972):
p(d; i(hc )) = 1 (hc)exp(-i(hc )d),
where d is the distance between ridges with height bigger than the cut-off height. i is the distribution's shape parameter i.e. the expected number of keels per kilometre. The expected distance between ridges is i• i is dependent on hc.
To obtain estimates of ship performance in ridged ice conditions, several random realizations of ridge field geometries are generated for each ice condition. Ship motion in the generated ridge fields is then simulated by the method described in Section 2.1 and the mean speeds are calculated for each ridge field. The probability of the ship becoming beset in ice is obtained by dividing the number of simulated cases where the ship gets stuck in ridges by the total number of cases. The output of the transit simulation is a distribution of mean
speeds and a probability of the ship becoming beset in given ice conditions.
Ridge field geometries are generated for the transit simulations by drawing random samples of distributions of ridge keel depths and spacings. Ridged ice fields are characterized by mean ridge keel depth, hk, which is used to parametrize Eq. (16), linear ridge density, l, used to parametrize Eq. (17) and level ice thickness, hi. The ridges are modelled as isosceles triangles with keel angle k. Ridges are generated until a pre-set length of ridge field is full. The first ridge starts at x = d, where d is a random distance distributed according to Eq. (17). If generated ridges overlap, the union of the overlapping keels is the used geometry; for completely overlapping ridges, the smaller ridge is thus discarded. A conceptual drawing of the ridge field model is shown in Fig. 3. If the ship encounters a situation where it does not have sufficient thrust to move forward or reverse, the simulation is marked as beset in ice and the mean speed is discarded from further analysis.
The sensitivity of the transit simulation method to the number of simulations per ice condition was tested. The mean of the mean speeds varied within 1.5% when the number of simulations varied between 100-1000. 200 simulations for each ice condition parametrized by {hi, hk,i} is used in this article.
3. Data-sets
Two data-sets were used to compare observed ship performance with simulation results. The first set is from the ice trials of S.A. Agul-has II and it is comprised of time-histories of ice thickness measured by an electromagnetic (EM) thickness sounding device and ship's speed. This data-set is described in Section 3.1.
The other data-set includes speed and location of a general cargo ship on a voyage between Kokkola and Vaasa in the Bay of Bothnia. The data are obtained from AIS transmissions of the ship. HELMI hindcast ice model as well as ice charts are used for information concerning ice conditions. This data-set is described in Section 3.2.
3.1. S.A.AgulhasII
Measurement results from the ice trials of S.A. Agulhas II were used for estimating of the resistance prediction method. The ice trials tookplace in the Bay of Bothnia on 21-22.3.2012.
S.A. Agulhas II is a polar supply and research vessel owned and operated by the South African Department of Environmental Affairs. She was built in Rauma shipyard by STX Finland in 2011-2012 and has an ice class of IACS PC5. The main characteristics of S.A. Agulhas II are presented in Table 4.
Agulhas was equipped with an EM device for measuring ice thickness for the ice trials. The method is described in detail in (Lensu et al., 2015). The device can not measure ice thickness exceeding 4 m reliably and it has a footprint of about 10 m. The difference between snow surface and ice bottom is obtained and observations on snow thickness are needed to calculate the ice thickness. (Lensu et al., 2015) Also, the power setting of the propulsion units, propeller pitch and revolutions, rudder angle, speed and course over ground and heading of the ship were measured. Furthermore, mechanical properties of ice were measured on both test days in the vessel's area of operations. (Bekker et al., 2014) The information on power setting as well as propeller pitch was used to distinguish a part of the data-set where the variations in speed are due to the ice conditions and not a result of crew action.
The flexural and compressive strength of ice were measured in the operating area of Agulhas II on three occasions. Measurements were carried out at 08:00 and 14:00 on 21.3.2012 and at 14:00 on 22.3.2012. The two measurements on 21.3. were conducted within about 6 km of each other and the last measurement was conducted about 100 km from the previous day's test sites. Crushing strength had a mean of 1.28 MPa with 0.38 MPa standard deviation. For flex-ural strength of ice, the mean and standard deviation were 0.4 MPa and 0.06 MPa respectively. (Suominen et al., 2013)
3.2. General cargo vessel
A data-set combining the speed and location of the vessel with relevant ice conditions has been created based on AIS-transmissions and the HELMI hindcast ice model (Montewka et al., 2015). Only parts of the voyage, when the ship was operating independently with full power, have been included in the data-set. For comparison with simulation results, two segments of the data-set were selected, during which the ice concentration was over 0.95 and the ship operated in the continuous mode i.e. did not have to ram ridges or become beset in ice.
The general characteristics of the ship are given in Table 5. The ship is designed to have the structure and power to navigate in difficult ice conditions unassisted. The draught of the ship at the time of the voyage was not known and the ship was assumed to operate in half-load conditions. This causes some uncertainty in the prediction of the performance of the vessel. (Montewka et al., 2015)
In the selected intervals of the voyage, the level ice thickness indicated by HELMI was 0.4 m. HELMI does not provide data on the small scale thickness variation of ridge fields but rather the total volume of ridged ice for larger areas (Haapala et al., 2005). Realizations of ridge field profiles were generated for the simulations using Eqs. (16) and 17. It is assumed that the mean sail height of ridges in the Baltic is constant at 0.6 m with a cut-off height of 0.4 m. (Lepparanta, 1981). With the assumption of constant mean ridge thickness, the equivalent thickness of ridged ice and ridge density are related by
heq = 0.022m2i ± 27%, (18)
Fig. 3. Conceptual drawing of the ridge field model. d is defined as the distance between adjacent ice thickness maxima. Sails and keels are assumed to coincide.
Table 4
Main characteristics of S.A. Agulhas II.
Ice class 1ACS PC5
DWT 5000 t
Length 135 m
Breadth 21.7 m
Draught 7.65 m
Power 9000 kW
Speed 14 kn
Table 5
Main characteristics of the bulk carrier used in the transit simulation.
1ce class IAS
DWT 21353 t
Length 149.3 m
Breadth 24.6 m
Draught 9.4 m
Power 9720 kW
Speed 15.5 kn
where heq is the equivalent thickness of ridged ice obtainable from the HELMI model. One ridge per km accounts on average for 2.2 cm of equivalent ice thickness. (Lepparanta, 1981) With the assumed triangular ridge model with Kk = 25° this corresponds to mean keel depth between 2.6-3.4 m. This means that the ratio of keel depth to sail height is 4.3-5.7. The range of the hk/hs ratio is somewhat larger than those given in Section 2.2 but well within reasonable bounds.
The intervals of the voyage selected for analysis had heq between 0.24-0.28 m. These values correspond to quite heavily ridged ice with ridge densities varying between 9-18 1/km. The ridge fields were generated for the transit simulations with the assumptions that one ridge per km corresponds to 16, 22 and 28 cm of equivalent ice thickness to compare the performance predictions over the range of the uncertainty in Eq. (18).
4. Sensitivity analysis
The simulation method contains many simplifying assumptions as well as parameters that are difficult or impossible to measure accurately in operational situations. Also, many properties of the ridged ice field are assumed to be constant but may vary significantly in nature. For these reasons, it is important to test the sensitivity of the model to variations in the parameter values. Patey and Riska (1999) performed a sensitivity analysis covering level ice thickness and ridging parameters. However, they did not check the sensitivity of the model to the underlying assumptions and had a different interpretation of the results. Further, they did not investigate the probability of a ship becoming beset in ice but rather selected simulated conditions to minimize the possibility of a ship becoming beset (Patey and Riska, 1999).
The excitation is composed of three resistance components and the net thrust. The resistance is the sum of level ice, Rl, resistance and resistance caused by displacing the ridge mass at the bow, Rb and sliding of the ridge mass along the midship, Rm.
The significance of the resistance components to the overall work was investigated by simulating 200 cases each of all of the combinations of hi e {0.30.40.50.6} m and heq e {0.050.10.150.20.25} m. The range of the equivalent ice thickness corresponds to roughly 2-11 ridges/km. Rl was the most significant component in all cases contributing 59-91% of overall work. Rb accounted for 8-38% and Rm for 0.5-2.5%. However, Rb can be the largest resistance component momentarily, accounting for up to 90% of the total resistance when the bow penetrates the deepest point of a ridge keel.
The sensitivity of the mode to the ice conditions is investigated by varying level ice thickness, hi, equivalent ice thickness, heq, as well as the coefficient of friction between ship and ice, ii. Also the thickness of the consolidated layer relative to the surrounding level ice thickness is varied as is the relation of heq and linear ridge density, i.
The assumptions of the model are tested by varying the coefficients Cp and Cm (Eqs. (5) and 7), that govern the magnitude of the ridge resistance at the bow and midship respectively. As shown in Eq. (6), the value of Cp is influenced by the porosity and internal friction angle of the ridge rubble and the density difference of the sea water and ice. Eq. (8) shows that Cm is influenced by the coefficient of friction between ship and ice, density difference of water and ice as well as Poisson's ratio and porosity of the ridge rubble. The uncertainty of Cp can be expressed as
d Cp d Pd
where the partial derivatives are evaluated at the values ofp, pA and 0r used in the baseline assumptions. The range of the uncertainty of the parameter values is set based on state-of-art knowledge.
(Mellor, 1980) The uncertainty is dominated by the internal friction angle of the ridge rubble.
Similarly, the uncertainty of Cp can be expressed as
The value of Cm is not dominated by a single parameter but the coefficient of friction has the greatest effect. The baseline assumptions, range of uncertainty of the parameters and relative significance of the parameters to the overall uncertainty is given in Table 6.
Also the distance that the ship reverses when ramming and the assumptions of the thrust model are included in the sensitivity analysis. There has been work done to determine the effect of ice covered water on propeller curves and on the propulsive coefficients of ships in ice. The studies deal with level ice and ice channels but not with situations, where the propeller is immersed in ice rubble (Kosti-lainen, 1981; Leiviska, 2004). To the knowledge of the authors, there are no data publicly available, which could be used to quantify the effect of ridge rubble to the propulsive efficiency of a ship. Also the difference in the propulsive efficiency when a ship operates forward or reverses cannot be accurately quantized with current knowledge. However, simulations are run with different assumptions to gain insight into the significance of thrust modelling on the simulation results.
The parameters and their values used in the sensitivity analysis are shown in Table 7, where the baseline case and all other test cases are listed. The baseline case corresponds to the assumptions made in Section 2 and the range of tested parameter values is decided based on literature (Hoyland, 2002; Lepparanta, 2005; Mellor, 1980).
All of the simulations of the sensitivity analysis use the general cargo ship described in Table 5. 200 simulations were conducted for each set of input parameter values to obtain a distribution of mean speeds and the probability of the ship becoming beset in ice.
4.1. Ice conditions
Level ice thickness has a clear effect on both the mean speed and probability of getting stuck in ice. Fig. 4 shows simulated mean speeds for variable level ice thickness with i =11 1/km, hr = 3 m and vi = 5 m/s. Expected values of mean speed decrease almost linearly with increasing ice thickness over the investigated range. The probability of getting stuck is also influenced significantly by the level ice thickness. When the ice thickness increases from 0.3 m to 0.6 m the probability of getting stuck in ice increases from 0.13 to 0.75.
Simulated mean speeds in varying ridge densities are shown in Fig. 5. hi = 0.4 m, hr = 3 m and initial speed is 5 m/s in all cases. There is a near linear decrease in the expected values of the mean speeds up to i =11 1/km but the trend is less clear for ridge densities of 11-13 ridges per km. The probability of becoming beset in ice increases from 0.06 to 0.47 within the range of investigated ridge densities.
Table 6
Breakdown of the uncertainty of the parameters Cp and Cm.
Relative effect
Parameter Baseline Uncertainty ACp ACm
0r 52.5° 11° 0.895
p 0.3 0.15 0.005 0.11
Pa 125 kg/m3 40 kg/m3 0.1 0.16
m 0.26 0.09 0.23
li 0.15 0.15 0.5
Table 7
The test matrix of the sensitivity analysis. hi [m] heq [m] Ch [m2] Cc i Cp [kg/m2s2]
0.4 0.25 0.022 1.5 0.15 7500
0.2 0.3 0.5 0.6 0.7
0.4 0.15
0.2 0.3
0.25 0.016
0.028 0.022
Cm [kg/m2 s2] t [s] Vi [m/s] Cr Cf Cth
45.9 900 5 1 1 1
1.3 1.8
1.5 0.05
0.1 0.2
0.15 3600 5700 7800 9900 12000
7500 35
45.5 56 66.5 77
45.9 1800
2700 3600
0.8 0.9
1 0.6 0.7 0.8 0.9
O— mean
---95% interval
.......min and max
0.4 0.5
h, [m]
2.5% qtl
8 10 ¡J, [1/km]
Fig. 4. Distributions of simulated mean speeds in varying level ice thickness. i =11 ridges/km and hr = 3 m. Initial speed 5 m/s. 200 simulations per ice thickness were performed to obtain the speed distributions.
Fig. 5. Distributions of simulated mean speeds in varying ridge densities. hi = 0.4 m and hr = 3 m. Initial speed 5 m/s. 200 simulations per ridge density were performed to obtain the speed distributions.
4.2. Modelling assumptions
The sensitivity of the model to the assumptions made in the modelling was tested by a sensitivity analysis. The varied parameters and their values are shown in Table 7. 200 simulations were performed for each test case. The same set of simulated ridge field
geometries were used for the sensitivity analysis concerning all of the parameters considered in this subsection.
A summary of the results of the sensitivity analysis is provided in Table 8. The most significant parameters with regard to model sensitivity are Cp, which determines the magnitude of the ridge resistance at bow and Ctf, which describes the decrease of propulsive efficiency
Table 8
Qualitative summary of the results of the sensitivity analysis.
Parameter Uncertainty Sensitivity mmean P(beset)
Cp High Medium High
Cm High Insignificant Insignificant
Cc Medium Insignificant Low
Ctf Medium Medium Medium
Ctb Medium Insignificant Insignificant
Cr High Insignificant Insignificant
li High Insignificant Insignificant
in ice conditions. Also, the initial velocity has a significant impact in the results.
The variation of Cp is very large when calculated using the whole range of state of art knowledge of internal friction angle of the ridge rubble. Consequently, the means of mean speeds vary from 15% higher to 8% lower compared to the baseline case when the value of Cp is varied between 3600-12,000 kg/m2 s2. The effect of the probability of becoming beset in ice is affected even more with the probability varying from 0 to 0.85. Fig. 6 shows the distributions of mean speeds and probability of becoming beset in ice over the range of Cp.
The ratio of net thrust in open water to the net thrust in ice conditions, Ctf, has a significant impact on both the mean speed and probability of becoming beset in ice. Ctf was varied from 0.7 to 1 and the mean of mean speeds was 22% lower for Ctf = 0.7 than for Ctf = 1. The probability of becoming beset in ice varied between 0.4 and 0.775. Fig. 7 shows the distributions of mean speeds and probability of becoming beset in ice over the range of Ctf.
The ratio of consolidated layer thickness to level ice thickness, Cc, does not have a significant effect on mean speed over the tested range of 1.3-1.8. The probability of becoming beset in ice is however affected. The probability of becoming beset in ice increases from 0.35 to 0.49 over the investigated range of Cc.
Also the initial speed of simulations and the length of the simulated time were varied to assess their effect on the behaviour of the model. Decreasing initial speed of the simulation from 5 m/s to
0.8 0.6 0.4 0.2 0
—e— Mean, v
-----95% qtl, V
...........Extrema, v
- X- - P(beset)
3600 5700 7800 9900 12000 cp [kg/m2s2]
Fig. 6. Distributions of mean speed and probability of becoming beset in ice for different values of Cp. Speeds are normalized to the mean of mean speeds of the baseline case.
1.2 1.1
0.9 0.8 0.7 0.6 0.5 0.4
—©— Mean, v
-----95% qtl, v
...........Extrema, v
- X - P(beset)
Fig. 7. Distributions of mean speed and probability of becoming beset in ice for different values of Ctf. Speeds are normalized to the mean of mean speeds of the baseline case.
2 m/s reduced the mean of mean speeds by 12%. The probability of becoming beset in ice increased from 0.4 to 0.55.
The length of the simulated time does not affect the mean of the mean speeds significantly but the variation in mean speeds is decreased with increasing simulation time. The probability of becoming beset in ice is increased with longer simulation times. With heavily ridged ice the effect is very large. With hi = 0.4 m/s, heq = 0.25 m and vi = 5 m/s, the probability of becoming beset in ice increased from 0.4 to 0.85 when the simulation time was increased from 900 to 3900 s.
5. Results
The model was tested for reasonable response to varying inputs and against full-scale measurements to determine its feasibility for the intended purpose.
The primary aim of the model is to provide reasonable estimate of ship speed and probability of getting beset in ice for ice conditions accounting for the presence of ice ridges.
Comparisons of simulations to full-scale data are presented in Section 5.1. Comparisons were made both for the resistance prediction method and the transit simulation model.
5.1. Comparison to full-scale data
5.1.1. S.A. Agulhas II
A time interval during which the ship operated in ridged ice with nearly constant power and minimal manoeuvring was selected for comparison with simulation results. The ice conditions were relatively easy for the selected interval. Rayleigh separation was used to identify ridge keels from the EM data to characterize the level of ridging. Ridges were determined to be local maxima of ice thickness separated by minima of depths no more than half of the thickness of the maximum. A cut-off depth of 1 m was used in the Rayleigh separation. Fig. 8 illustrates the process of Rayleigh separation.
Over the 3.1 km of ridged ice there were 13 ridges of depth greater than 1 m, leading to a ridge density of 4.1 ridges per km with mean depths of 2.25 m. There was a 500 m stretch of ice nearby with mean thickness of 0.3 m in the data and this was assumed to be the
<D 1.5 c
50 100 150 200 250 300 350 x[m]
..............Cut-off thickness
O Ridge
x Non-Rayleigh maxima
■ Below cut-off maxima
4 2 0 -2
v . [m/s]
sim 1 '
v [m/sl
meas 1 J
-z[m] -© ridge
0 500 1000 1500 2000 2500 3000 x[m]
Fig. 8. Rayleigh separation illustrated. All local maxima of the thickness profile are highlighted. Ridges fulfil the Rayleigh criteria, non-Rayleigh maxima do not have surrounding minima of more than half their height and below cut-off maxima are below the set cut-off thickness.
Fig. 9. Simulated and observed time histories of the speed of S.A. Agulhas II operating in ridged ice. Circles indicate the deepest points of Rayleigh-separated ridges with cut-off depth of 1 m. z is the bottom profile of the ice field measured by an EM device.
level ice thickness for the whole area. The thickness of consolidated layers was set to 0.45 m based on the assumption that the thickness of consolidated layers is 1.5 times the thickness of the surrounding level ice.
The strength of ice was set for the simulations based on the measurements conducted during the ice tests described in Section 3.1. Table 9 shows the parameter values used for the simulations. Initial velocity for the simulation was the same as that in the comparison data-set.
Fig. 9 shows the simulated and observed speed of Agulhas and the ice thickness profile during the studied interval. The ridges defined by Rayleigh separation are also indicated in the figure. There is a very good correspondence between the simulated and observed mean speeds. The simulated mean speed is 2.5% higher than the observed and the speed profiles show qualitatively similar behaviour. The simulated speed varies within 3.8-4.8 m/s while the observed speed has a slightly larger variation of 3.1-5.2 m/s.
5.1.2. Transit simulation
Two time intervals of the voyage of the general cargo carrier, 12:45-13:44 and 14:40-15:13, were selected for the comparison with simulated results. Criteria for the data were total ice concentration of over 95% indicated by HELMI and independent navigation. The ship operated in difficult ice conditions during the selected time intervals.
HELMI model results were used to determine the ice conditions encountered by the ship. Level ice thickness was rounded to the nearest 10 cm. There were four different segments in three distinct ice conditions. Level ice thickness was 0.4 m in all of the cases. During 12:45-13:44, there were three different segments with equivalent
Table 9
Mechanical properties of ice used for the simulation of S.A. Agulhas II in ridged ice.
sf 0.3 MPa
s c 1.4 MPa
E 5 GPa
li 0.15
h¡ 0.3 m
Pw 1025 kg/m3
Pi 900 kg/m3
ice thickness of 0.24 m, 0.26 m and 0.28 m. Between 14:40-15:13, the equivalent ice thickness was 0.28 m. Three sets of 200 simulations were run for each segment corresponding to the average and extreme assumptions of Eq. (18). The evaluated ridge densities, rounded to the nearest integer value, were 9-15 1/km, 9-17 1/km and 10-18 1/km for equivalent ice thickness 0.24 m, 0.26 m and 0.28 m respectively.
Fig. 10 shows the observed mean speeds and distributions of simulated mean speeds for the four investigated segments. Simulations have the same initial speeds as the time histories used to determine the observed average speeds. The observed mean speeds are mostly reasonably close to the expected values of the simulations. For the fourth case, with heq = 0.286 m and vi = 3.6 m/s the observed mean speed is notably higher than the simulations would suggest. This is also the case with the highest observed mean speed, even though the ice conditions are the most severe. This can be a consequence of the natural variation of mean speeds due to small-scale variation in the ice conditions or it may be due to the presence of leads or channels not included in the ice data.
The different assumptions in mapping equivalent ice thickness to ridge density have a significant impact on the predicted mean speeds. Also, the predicted probabilities of the ship becoming beset in ice are significantly affected. The probabilities varied between
Observed mean speed
h =0.242 m h =0.264 m h =0.286 m h =0.286 m eq eq eq eq
¥¡=3.6 m/s ¥¡=4.8 m/s Vj=2.8 m/s ¥¡=3.6 m/s
Fig. 10. Boxplot of simulated mean speeds of the general cargo vessel in ridged ice with observed mean speeds of four legs of the journey. Three sets of 200 simulations have been calculated for each case, corresponding to different assumptions of ridge density. hi = 0.4 m in all cases. The bull's-eye indicates the mean value of simulated speeds, the thicker line covers the 50% central interval and the narrow line covers the range of the data except for outliers, which are indicated by small circles.
0.17-0.61, 0.26-0.67, 0.43-0.81 and 0.35-0.78 for the four cases. Significant probabilities of getting stuck in ice are in line with the ship becoming beset twice in similar ice conditions during the journey. The results are summarized with the corresponding ice conditions in Table 10.
6. Discussion
The sensitivity analysis showed that parameter values have to be very carefully chosen for certain parameters. The model is especially sensitive to Cp, which is dominated by the internal friction angle of the ridge rubble. This is problematic as the friction angle is difficult to measure in full scale and the accepted values for it cover a large interval.
Also, the assumptions made in the modelling of propeller thrust in ridge rubble were shown to have a significant effect on the model predictions. Especially the effect of ridge rubble in the propulsive efficiency is significant while the effect of reversing on the propulsive efficiency has no significant impact on the results.
Furthermore, the thickness of the consolidated layers does not affect the speed predictions but has a noticeable effect on the probability of becoming beset in ice.
The sensitivity analysis also showed that the length of the simulated transit in ridges has a significant effect on the probability of the ship becoming beset in ice and in the variance of the predicted mean speeds. Also,the initial velocity affects both the mean speeds and probability of becoming beset. When the method is used for transit simulations, care has to be taken to determine properly the length of route segments with similar ice conditions and to use the speed predictions of previous segments as a basis for setting the initial speeds.
Ridged ice is typically parametrized with expected values of ridge sail heights and spacings. Parametrization can also be made for ridge keels but in operational situations there are typically no direct observations of ridge keels. Inferences on ridge keel sizes and spacings can be made based on observed relations of ridge sail and keel occurrence and sizes.
For ridge fields parametrized with i, hr and ht, the model predicts a large scatter of mean speeds. The limited observations available support the idea that performance in given ice conditions should be described as a distribution of speeds rather than a point estimate but it is not possible to deduce from the data whether the simulated variation is too large. The classical models of ridge occurrence assume that ridge spacings are independently distributed. Lensu (2003) proposed a model in which ridge spacings are correlated and ridges are more likely to be situated near other ridges. A part of this large variability could be explained by clustering effects. The current method of generating ridge field profiles does not model clustering explicitly but realizations of ridge fields with varying levels of clustering are produced. Possible correlations between measures of clustering and predicted mean speeds could give an indication of the importance of ridge clusters for the model.
Another significant source of uncertainty is the mapping from equivalent ice thickness to ridge density. In the current form, this causes very significant differences in the performance predictions.
Table 10
Summary of the transit simulation results and corresponding ice conditions. h, = 0.4 m in all cases, hk varies between 0.5 m and 2.8 m with a mean of 1.5 m.
l [1/km] Vi [m/s] Vm [m/s] Vobs [m/s] P(beset)
9-15 3.6 3.3-3.9 3.2 0.17-0.61
9-17 4.8 3.5-4.0 3.2 0.26-0.67
10-18 2.8 3.0-3.7 3.5 0.43-0.80
10-18 3.6 3.1-3.7 4.3 0.35-0.78
The parameters for the distributions of ridge sizes and spacings are scale dependent, i.e. they vary depending on how large areas are observed in order to determine the parameters. The correct scale used for determining the parameter values is not readily apparent and should be studied further. The parameters are also dependent on the cut-off value for sail height or keel depth that is used when the parameters are determined from data. (Lepparanta, 2005).
Also, modelling the amount and spatial distribution of underwater ice rubble based on the characteristics of the ice cover observable from the surface is somewhat uncertain. For large length scales, the correspondence is clearer but the importance of small scale phenomena is not clear at the moment. (Wadhams, 1980)
Modelling of ridge resistance in this model and the superposition of the different resistance components contains significant simplifying assumptions. No speed dependence is modelled for ridge resistance and the model is meant for low speeds (Malmberg, 1983). However, the results of comparing observed performance of S.A. Agulhas II and the simulation results is encouraging, although data is only available for quite easy ice conditions.
Future improvements of the model should include at least an improved ridging model, that takes ridge clustering explicitly into account. Also, the mapping from information available in ice forecasts and remote observations to the parameters of ridge field model used for the simulations should be further investigated. The current thrust model is very simple and more realistic modelling of the ship's thrust has a potential for improving the reliability of model results. Also, the behaviour of the ship in ramming ice ridges should be better modelled. Finally, the effect of compressive ice is important for navigation and including its effects would greatly increase the predictive power of the model.
7. Summary
A method for estimating a ship's resistance in ridged ice was reviewed and used as a basis for a novel method for predicting the performance of ships in ridged ice conditions. Performance is understood here to mean attainable average speed and ability to operate independently. Ice condition is parametrized by level ice thickness, ridge density and mean ridge size.
Multiple deterministic simulations of a ship transiting ridged ice are performed for randomly generated ridge-field profiles to obtain a distribution of mean speeds for an ice condition. The approach captures some key elements of navigating in ridged ice that have not been visible in previous modelling schemes. The considerable variability of mean speeds within an ice condition is reflected also in the available data and estimations of probability of the ship becoming beset in ice have not been included in previous physically based methods.
A sensitivity analysis found that the model is sensitive to assumptions concerning the physical properties of the ridge rubble and thrust in ice rubble. Also the importance of using simulations of suitable length and with realistic initial velocities was shown.
The suitability of the resistance models for predicting ship's speed was shown for lightly ridged conditions by comparing observed and simulated time histories of ship speed. Ability of the model to make meaningful predictions for different ice conditions was also shown by comparisons with full-scale data. Although the used data-sets are not very large, the results are good.
The presented approach attempts to advance the ship in ice transit modelling, by improving the ridge field model. The obtained results demonstrate that the proposed model delivers reasonable results in terms of ship speed and the probability of getting stuck. This makes it possible to apply the model for practical use for instance for route planning in ice or estimating the operability of a ship in given ice conditions.
Acknowledgements
This work has received funding from the European Commission (605190) funded project "Joint Operation for Ultra Low Emission of Shipping" (JOULES) as well as the academy of Finland projects "Kara Arctic Monitoring and Operation Planning Platform" (KAMON) and "Variation of Antarctic sea ice thickness and its effect on the load level of ice navigating" (ANTLOAD) as well as TEKES project "Vessel Operations and Routing in Ice Conditions" (VORIC).
Justin Champion read the manuscript and gave valuable comments on English grammar and usage.
References
Abdelnour, R., Comfort, G., Peirce, T., 1991. Single pass ridge penetration model. Proceedings of the 11th International Conference on Port and Ocean Engineering under Arctic Conditions. vol. 2. pp. 600-622.
Bailey, E., Taylor, R., Croasdale, K., 2015. Mechanics of ice rubble over multiple scales. OMAE2015.
Bekker, A., et al. 2014. Full-scale measurements on a polar supply and research vessel during maneuver tests in an ice field in the Baltic sea. Proceedings of the 33rd International Conference on Ocean, Offshore and Arctic Engineering.
Bowen, R., Topham, D., 1995. A study of the morphology of a discontinuous section of a first year arctic pressure ridge. Cold Reg. Sci. Technol. 24 (1), 83-100.
Guinnes, R., et al. 2014. A method for ice-aware route optimization. Position, Location and Navigation Symposium-PLANS 2014. pp. 1371-1378.
Haapala, J., Lonnroth, N., Stossel, A., 2005. A numerical study of open water formation in sea ice. J. Geophys. Res. 110, C09011.
Hibler, W., Weeks, W., Mock, S., 1972. Statistical aspects of sea-ice ridge distributions. J. Geophys. Res. 77 (30), 5945-5970.
Hopkins, M., Hibler, W., Flato, G., 1991. On the numerical simulation of the sea ice ridging process. J. Geophys. Res. 96 (C3), 4809-4820.
Hoyland, K., 2002. Simulations of the consolidation process in first-year sea ice ridges. Cold Reg. Sci. Technol. 34 (3), 143-158.
Kamesaki, K., Kishi, S., Yamauchi, Y., 1999. Simulation of NSR shipping based on year-round and seasonal operation scenarios. INSROP Working Paper 164-1999. INSROP.
Kankaanpaa, P., 1989. Structure of first year pressure ridges in the Baltic sea. Finnish Institute of Marine Research.
Kankaanpaa, P., 1997. Distribution, morphology and structure of sea ice pressure ridges in the Baltic sea. Fennia 175 (2), 139-240.
Keinonen, A., 1979. An Analytical Method for Calculating the Pure Ridge Resistance Encountered by Ships in First Year Ridges. Helsinki University of Technology. (Ph.D. thesis).
Kostilainen, V., 1981. Performance of marine propellers in ice-clogged channels. Tech. Rep.. Winter Navigation Research Board.
Kotovirta, V., et al. 2009. A system for route optimization in ice-covered waters. Cold Reg. Sci. Technol. 55, 1.
La Prairie, D., et al. 1995. A transitsimulation model for ships in Baltic ice conditions. Tech. Rep.. Helsinki University of Technology, Ship Laboratory.
Leiviska, T., 2004. Laivan propulsiokertoimet jaissa ja avovedessa (propulsive coefficients of a ship in ice and open water). Tech. Rep. M-287. Helsinki University of Technology.
Lensu, M., 2003. The Evolution of Ridged Ice Fields. Helsinki University of Technology, Ship Laboratory, Espoo, Finland. (Ph.D. thesis).
Lensu, M., et al. 2015. Measurements of antarctic sea ice thickness during the ice transit of S.A. Agulhas II. Proceedings of the 23 rd International Conference on Port and Ocean Engineering under Arctic Conditions.
Lepparanta, M., 1981. On the structure and mechanics of pack ice in the Bothnian Bay. Tech. rep. 248. Finnish Marine Research.
Lepparanta, M., 2005. The Drift of Sea Ice. Springer.
Lepparanta, M., et al. 1995. The life story of a first-year sea ice ridge. Cold Reg. Sci. Technol. 23(3), 279-290.
Lewis, J., Lepparanta, M., Granberg, H., 1993. Statistical properties of sea ice surface topography in the Baltic Sea. Tellus 45 A, 127-142.
Malmberg, S., 1983. Om fartygs fastkling i is (of ship's becoming beset in ice). Helsinki University of Technology. In Swedish MA thesis.
Melling, H., Topham, D., Riedel, D., 1993. Topography of the upper and lower surfaces of 10 hectares of deformed sea ice. Cold Reg. Sci. Technol. 21, 349-369.
Mellor, M., 1980. Ship resistance in thick brash ice. Cold Reg. Sci. Technol. 3, 305-321.
Montewka, J., Goerlandt, F., Kujala, P., 2015. Towards probabilistic models for the prediction of a ship performance in dynamic ice. Cold Reg. Sci. Technol. 112, 14-28.
Mulerhin, N., Eppler, D., Sodhi, D., 1999. A time and cost prediction model for northern sea route shipping. INSROP Working Paper 150-1999. INSROP.
Newmark, J., 1959. A method of computation for structural dynamics. Proceedings of the American Ociety of Civil Engineers vol. 85 (3), 67-94.
Parmerter, R., Coon, M., 1972. Model of pressure ridge formation in sea ice. J. Geophys. Res. 77 (33), 6565-6575.
Patey, M., Riska, K., 1999. Simulation of ship transit through ice. INSROP, INSROP Working Paper 155-1999.
Riska, K., et al. 1997. Performance of merchant vessels in ice in the Baltic. Tech. Rep. Report 52. Winter Navigation Research Board, Helsinki.
Strub-Klein, L., Sudom, D., 2012. A comprehensive analysis of the morphology of first-year sea ice ridges. Cold Reg. Sci. Technol. 82, 94-109.
Suominen, M., et al. 2013. Full-scale measurements on board PSRV S.A. Agulhas II in the Baltic sea. Port and Ocean Engineering under Arctic Conditions.
Timco, G., Burden, R., 1997. An analysis of the shapes of sea ice ridges. Cold Reg. Sci. Technol. 25(1), 65-77.
Tin, T., Jeffries, M., 2003. Morphology of deformed first-year sea ice features in the southern ocean. Cold Reg. Sci. Technol. 36 (1-3), 141-163.
Tuhkuri, J., Lensu, M., Saarinen, S., 1999. Laboratory and field studies on the mechanics of ice ridge formation. Proc. POAC 1999.
Veitch, B., et al. 1991. Field observations of ridges in the northern Baltic sea. Proc. POAC 1991.
Wadhams, P., 1980. A comparison of sonar and laser profiles along corresponding tracks in the arctic ocean. Sea Ice Processes and Models. pp. 283-299.
Wadhams, P., Davy, T., 1986. On the spacing and draft distributions for pressure ridge keels. J. Geophys. Res. 91 (C9). pp. 10,697-10,708.
Wadhams, P., Horne, R., 1980. An analysis of ice profiles obtained by submarine sonar in the Beaufort sea. J. Glaciol. 25 (93), 401-424.
Wright, B., Hnatiuk, J., Kovacs, A., 1978. Sea ice pressure ridges in the Beaufort sea. Proceedings of IAHR Ice Symposium.
