Maxime Bavencoffe
email: maxime.bavencoffe@univ-tours.fr

Marc Lethiecq
email: marc.lethiecq@univ-tours.fr

A Loukkal

M Lematre

Modeling and numerical study of the influence of imperfect interface properties on the reflectance function for anisotropic multilayered structures

Keywords: anisotropic multilayer, imperfect interface, reflectance function, guided waves, V(z) method

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The production of adhesively bonded structures has significantly increased in several industries such as aerospace, automotive or marine. The great interest for adhesive bonding is due to its low cost, low weight, the good distribution of stress it insures and 2 the ease of assembly. Nevertheless, adhesive bond defects can appear. They can be mainly classified in two different types: adhesive defects related to the weak bond between the adhesive and the adherent, and cohesive defects related to the mechanical properties of the adhesive [START_REF] Santos | Ultrasonic guided waves scattering effects from defects in adhesively bonded lap joints using pitch and catch and pulse-echo techniques[END_REF]. The adhesion degradation can be caused or accelerated by moisture infiltration, through the adhesive layer or along the interface, and improper assembly process due to contamination and micro-cracks introduced at the interfaces. The ability of ultrasound waves to penetrate opaque media, with possibility to generate guided waves, makes them the most promising in order to assess the adhesion quality between coatings and substrates, two substrates or two layers. Many authors have studied the interaction between ultrasonic waves and adhesive layers [START_REF] Drinkwater | A study of the interaction between ultrasound and a partially contacting solid-solid interface[END_REF][START_REF] Rokhlin | Ultrasonic characterization of plasma spray coating[END_REF][START_REF] Guo | Modeling and acoustic microscopy measurements for evaluation of the adhesion between a film and a substrate[END_REF][START_REF] Baltazar | Inverse ultrasonic determination of imperfect interfaces and bulk properties of a layer between two solids[END_REF][START_REF] Siryabe | Apparent anisotropy of adhesive bonds with weak adhesion and non-destructive evaluation of interfacial properties[END_REF]. Their works addressed the case of normal and oblique incidence waves in pulse echo mode [START_REF] Drinkwater | A study of the interaction between ultrasound and a partially contacting solid-solid interface[END_REF][START_REF] Baltazar | Inverse ultrasonic determination of imperfect interfaces and bulk properties of a layer between two solids[END_REF] and transmission mode [START_REF] Siryabe | Apparent anisotropy of adhesive bonds with weak adhesion and non-destructive evaluation of interfacial properties[END_REF]. They evaluated also the adherence quality using Rayleigh waves [START_REF] Rokhlin | Ultrasonic characterization of plasma spray coating[END_REF] and leaky SAWs [START_REF] Guo | Modeling and acoustic microscopy measurements for evaluation of the adhesion between a film and a substrate[END_REF].

In this work, the interest is directed towards multilayered structures. The ultrasonic characterization of interfaces in this specific case has been assessed by several researchers [START_REF] Siryabe | Apparent anisotropy of adhesive bonds with weak adhesion and non-destructive evaluation of interfacial properties[END_REF][START_REF] Wang | Evaluation of interfacial properties in adhesive joints of aluminum alloys using angle-beam ultrasonic spectroscopy[END_REF][START_REF] Vlasie | Mechanical and acoustical study of a structural bond: comparison theory/numerical simulations/experiment[END_REF][START_REF] Ismaili | Determination of epoxy film parameters in a three-layer metal/adhesive/metal structure[END_REF][START_REF] Budzik | Monitoring of crosslinking of a DGEBA-PAMAM adhesive in composite/aluminium bonded joint using mechanical and ultra-sound techniques[END_REF]. For example, Vlasie et al. [START_REF] Vlasie | Mechanical and acoustical study of a structural bond: comparison theory/numerical simulations/experiment[END_REF] proposed an interface model simulating the damage evolution of a structural joint. The originality of their model lies in the fact that it allows to access the stiffnesses of the interface that cannot be derived from mechanical tests. To do so, they modeled the adhesive layer by a geometrical interface with a uniform mass/springs distribution. Then, they performed a comparison between the cutoff-frequencies of the ultrasonic guided modes for a trilayer model and a mass/springs model of the adhesive layer. This comparison allows to derive the longitudinal and transversal stiffnesses with respect to frequency. Thus, a comparison between simulation and experimental measurements enabled them to validate their model and evaluate the elastic parameters. Ismaili et al. [START_REF] Ismaili | Determination of epoxy film parameters in a three-layer metal/adhesive/metal structure[END_REF] proposed a method allowing to determine the thickness, the longitudinal and shear waves velocities of an adhesive layer embedded in a 3-layer structure. They used, subsequently, these parameters in order to study the sensitivity of some guided modes with respect to the parameters in question.

In this paper, we will study the influence of different parameters of an adhesive interface layer on the resonance modes (minima of the modulus of the reflectance function) of a multilayered structure. This study is included in the frame of V(z) measurements in order to characterize adherence properties in multilayered structures. When inversing the V(z) experimental data, it allows to rebuild the reflectance function R(θ) with a simple experimental set-up, that works in normal incidence only, with the use of a forward transducer of large angular aperture. The guided modes of the structure are sensitive to the interface conditions between layers, for this, each guided mode needs to be isolated. Thus, the reflectance function is used in this study as a parameter of control since it gives the resonance modes which are shown to correspond to the guided modes of the structure, and thus, it avoids the difficulty in guided waves measurements, thanks to the V(z) method as previously described.

An imperfect interface model inspired by Rokhlin's imperfect anisotropic interface layer [START_REF] Rokhlin | Ultrasonic wave interaction with a thin anisotropic layer between two anisotropic solids: Exact and asymptotic-boundary-condition methods[END_REF] has been implemented to model the reflectance function of an anisotropic multilayer possessing an interface layer offering a variable adherence quality. Rokhlin simulates a fictitious porous interface layer, created from a Hill's approximation of the substrate underneath the interface, with the possibility of changing the porosity factor and the orientation of the pores. In our case, the solid matrix of the porous medium is made of an epoxy resin and the adhesive layer thus obtained is embedded in a multilayered structure and not between two semi-infinite media such as in the case addressed by Rokhlin. From a practical point of view, this model allows for the simulation of an altered interface possessing defects such as delamination and porosity accounted for by the vacuum inside the pores. The possibility to vary the orientation of the pores allows also to investigate such defects characterized by a privileged orientation within the adhesive layer.

In this paper, the theoretical background, based on the transfer matrix method [START_REF] Brekhovskikh | Waves in Layered Media[END_REF] and the porous interface layer, is described first. Afterwards, we establish the equivalence between the resonance modes obtained for a bilayer (containing an interface layer) immersed in water and the dispersion curves of the guided waves calculated for the same structure in vacuum. In the same section, the waves coupling phenomenon between the different layers is highlighted. Finally, the influence of different properties, such as thickness and porosity, of the adhesive layer on the resonance modes of the multilayer is studied as a function of the frequency and the incident angle.

Theoretical background

Global transfer matrix of a multilayer

An anisotropic multilayered structure consisting of N layers is considered (Figure 1). Each layer k of this structure has a thickness d and the total thickness is equal to D. The multilayer structure is immersed in water (media 0 and N + 1) and a longitudinal wave is assumed to be incident on the structure with an angle θ . Due to the anisotropic nature of the materials and using the superposition principle, in each layer, six waves can propagate in the medium [START_REF] Hosten | Transfer matrix of multilayered absorbing and anisotropic media. Measurements and simulations of ultrasonic wave propagation through composite materials[END_REF]. Indeed, in an anisotropic medium, a quasilongitudinal and two quasi-shear waves can propagate, and the fact that these waves are reflected by the bottom surface of the layer induces that these three kinds of waves propagate also in the opposite direction of the incident ones. The incident wave is reflected on the multilayer with the same angle θ , previously introduced, with respect to the normal as shown on Figure 1. Using this feature, we can calculate the reflectance function related to the multilayered structure.

Following the well-known Thomson-Haskell matrix method [START_REF] Lowe | Matrix techniques for modeling ultrasonic waves in multilayered media[END_REF][START_REF] Rose | Ultrasonic waves in solid media cambridge university press[END_REF], displacements and stresses existing at the bottom of a layer are connected to those at the top of the same layer by means of a layer transfer matrix as described in Eq. (1).

= ′ ′ (1)
where , , and , , are respectively the components of the displacements of the lower and the upper face of a layer. Similarly, , σ , σ and σ , σ , σ are the stresses applied on the lower and upper faces of a layer, respectively. is the 6*6 layer transfer matrix [START_REF] Hosten | Transfer matrix of multilayered absorbing and anisotropic media. Measurements and simulations of ultrasonic wave propagation through composite materials[END_REF]. The continuity conditions then make it possible to establish a global transfer matrix resulting from the product of the transfer matrices of the successive layers. This global matrix allows the displacements and stresses of the top layer to be related to those of the bottom one as follows:

= ! " # $ (2)
where is the global transfer matrix of the multilayered structure, is the interface matrix relating the displacements and stresses of two consecutive layers (k) and (k+1), is the local transfer matrix of layer k and $ is the interface matrix relating the displacements and stresses of the upper layer of the multilayer to those existing in the coupling fluid (water in our case).

The interface between two layers has different properties from those of the effective layers forming a multilayered structure. The nature of the interface depends on the type of solid contact made (bonding, various metallurgical processes, etc.). Conventional interface conditions that replicate a perfect and infinitely thin interface are not suitable for describing real interfaces, which are often not perfect. Hence, an imperfect anisotropic interface layer model is introduced here. This imperfect interface is involved in the establishment of the global transfer matrix of the multilayered structure via the interface matrices named that were introduced in Eq. [START_REF] Drinkwater | A study of the interaction between ultrasound and a partially contacting solid-solid interface[END_REF]. The interface layer used in this work does not take into account specifically cohesive or adhesive defects but describes globally the strength of a bond between two layers. This is due to its small thickness in comparison with the thicknesses of the metallic layers and the ultrasonic wavelength that allows us to consider it as a very thin interface layer between two media.

Modeling of the interface layer

When an interface is imperfect and has different micro-defects, it can be considered as a multi-phase layer characterized by effective elastic properties of its own [START_REF] Rokhlin | Analysis of ultrasonic wave interaction with imperfect interface between solids[END_REF]. The imperfect interface is modeled here by a thin orthotropic layer containing cylindrical pores parallel to the interface [START_REF] Rokhlin | Ultrasonic wave interaction with a thin anisotropic layer between two anisotropic solids: Exact and asymptotic-boundary-condition methods[END_REF]. This virtual layer consists of an isotropic epoxy matrix in which pores (considered as vacuum) are integrated. The effective elastic constants of this layer are calculated using Christensen's two-phase composite model [START_REF] Christensen | Mechanics of Composite Materials[END_REF]. The pores direction can be shifted by an azimuthal angle φ with respect to the incident plane. The five elastic constants & , & , & , & and & '' refer to the five effective independent constants that fully describe the stiffness tensor of the interfacial medium. Christensen's model aims to provide analytical expressions for these five constants by taking into account the properties of the pores, the matrix and the respective volumes occupied by each of these phases. To do this, a geometric model of the composite material must be introduced. The most commonly used model is the cylinder composite model developed by Hashin and Rosen [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF]. The fibers (which correspond to the empty pores in our case) are made of infinitely long circular cylinders incorporated in a continuous matrix, as described in Figure 2. [START_REF] Christensen | Mechanics of Composite Materials[END_REF] In order to obtain the effective elastic constants of the interface layer, the expressions of the moduli ( , ) , * , + and + are required. ( is the effective uniaxial Young's modulus, ) is the effective Poisson's coefficient where x1 axis is in the direction of the fiber (pore), * is the bulk modulus, + is the shear modulus in the fiber direction and finally + is the transverse shear modulus in the direction perpendicular to that of the pores, i.e. in the x direction. The first four parameters are expressed as follows: )A> (3.a) where P denotes the porosity factor defined by the ratio of the volume occupied by the pores to that occupied by the solid material, ( . and ( 2 are the Young's moduli of the materials constituting respectively the fibers (vacuum in our case) and the matrix (epoxy resin in our case), ) . and ) 2 are the respective Poisson's coefficients of the materials constituting the fibers and the matrix, F . and F 2 are the respective bulk moduli of the materials constituting the fibers and the matrix, and finally + . and + 2 are the respective shear moduli of the materials constituting the fibers and the matrix. The elastic parameters of the epoxy matrix are calculated using the data available in a paper published by Rokhlin et al. [START_REF] Rokhlin | Ultrasonic evaluation of environmental durability of adhesive joints[END_REF], as for the parameters related to the vacuum filling the pores, we set all of these parameters to an infinitely small value (eps function in Matlab ® ).

( = -( . / (1 1 -)( 2 /
) = (1 1 -)) 2 / -) . /
In order to obtain the value of + the second degree equation below must be solved:

M I : 9@ : 8 J / 2O I : 9@ : 8 J / & = 0 (4)
& = + / * (5.c) & = 1+ / * (5.d) & '' = + (5.e)
In this model, the alteration of the quality of adherence is directly related to the increase of the porosity in the interface layer. Since the quality of an interface is dependent on the shear modulus (whose value is responsible for the transmission of the shear components of stresses and displacements), it is crucial, in our case, to know the evolution of the shear moduli + and + with respect to the porosity factor P. Figure 3 shows the evolution of both shear moduli, defining the interface layer, with respect to the porosity. The value of these moduli decreases as the porosity grows until, logically, it reaches 0 when P = 1. The porous interface layer described here is included in the global transfer matrix method through a 6*6 matrix. The elements of this interface matrix are obtained in the same way as the local transfer matrix calculated for the metallic layers. Indeed, the difference compared to the other local transfer matrices of the structure lies in the thickness of the interface layer, which is much lower than the other thicknesses involved. The calculation of the velocities of the waves propagating in the medium and their polarizations are performed using the elastic constants expressed in Eq. ( 5), as well as the density U V of the interface material. The latter is obtained by the relation:

U V = (1 1 -)U WXY (6)
Where U WXY is the density of the epoxy resin, fully polymerized, acting as a solid matrix for the fictitious porous interfacial medium.

Free guided waves and reflectance function of the multilayer

When writing the continuity conditions of displacements and stresses at the top and bottom interfaces with those of the surrounding fluid and using Eq. (1) and Eq. ( 2), one can obtain the expression of the reflection coefficient as follows: [START_REF] Wang | Evaluation of interfacial properties in adhesive joints of aluminum alloys using angle-beam ultrasonic spectroscopy[END_REF] where c is the angular frequency, d is the acoustic impedance of the coupling fluid and the e V` are the determinants of matrices made of combinations of coefficients of the global transfer matrix . The expressions of e V` are given in Appendix B.

Z = [ 9G (\Y]^_) 9 5`ab _ \Y]^_[ 99 >ab _ (`[ GG \Y]^_>ab _ [ G9 ) [ 9G (\Y]^_) 9 >`ab _ \Y]^_[ 99 >ab _ (`[ GG \Y]^_5ab _ [ G9 )
In order to analyze the wave propagation inside the multilayered structures, the free guided waves dispersion curves have to be calculated. If we consider a multilayered structure surrounded by vacuum, it is possible to calculate the dispersion curves of the guided waves by adopting the free surfaces conditions. These boundary conditions result in the following relation between the displacements and stresses located respectively at the bottom (referred to by the superscript bot) and the top (referred to by the superscript top) of the structure: As we want to determine the free guided waves, corresponding to the multilayered structure surrounded by vacuum, it implies that the stress components , and at the top and bottom interface of the multilayer have to be null.

Taking into account the nullity of the previous stress components, Eq. ( 8) possesses non trivial solutions only if:

p m 3 m 3 m 3 m n m n m n m ' m ' m ' p = 0 (9)
The numerical resolution of Eq. ( 9) thus allows the phase velocities of the free guided waves to be obtained with respect to frequency. This equation can also be established by nullifying the denominator of the reflection coefficient, presented in Eq. ( 7), when setting the parameter d to zero. Indeed, since the poles of the reflection coefficient correspond to the generalized guided waves radiating in the incident fluid environment ("leaky guided waves"), the free guided waves are obtained by setting the acoustic impedance d of the incident coupling fluid to zero, leading to equivalently determine the roots of the scalar value e . Thus, when looking at the expression of e in Appendix B, the equality e = 0 strictly corresponds to the equation [START_REF] Ismaili | Determination of epoxy film parameters in a three-layer metal/adhesive/metal structure[END_REF]. The solutions of Eq. ( 9) are obtained numerically using the MATLAB ® software. In our case, for each frequency, the phase velocities of the guided waves are derived via the process described above. These results are plotted, with respect to frequency, in order to build the dispersion curves of the free guided waves. This protocol provides accurate dispersion curves, even if some modes appear discontinuous due to the computation steps of the frequency and incident angle.

Reflectance function minima and free guided waves dispersion curves for a perfect interface

In Section 2, the methodology to obtain the free guided waves and the reflectance function has been exposed.

The multilayered structure that will serve as a reference and point of comparison for our simulations is shown in Figure 4, and consists in an assembly of a 1.2 mm thick cubic anisotropic nickel (Table 1) layer and a 0.5 mm thick cubic anisotropic copper (Table 1) layer connected by a 2 μm thick interface layer. The interface layer, here, is considered as perfect (no porosity) and made of an epoxy resin (Table 1).

Validity of the guided waves dispersion curves and their comparison with the reflectance function

The multilayer is considered infinite along the x and x axes and its surfaces are free (structure placed in vacuum). The dispersion curves of the free guided waves can be plotted in an equivalent way as a function of phase velocity or incidence angle. Indeed, by identifying them to the corresponding leaky guided modes, and according to Snell's law, the angle of incidence { corresponding to each free guided mode is given by:

{ = arcsin • ' _ \ ƒ " ( 10
)
where … is the wave velocity in water and † ‡ is the phase velocity of the guided wave.

The representation of { as a function of frequency facilitates comparison between the reflectance function amplitude and the dispersion curves of the free guided waves. Thus, it should be noted that a higher phase velocity † ‡ implies a smaller incidence angle { , at a given frequency. The dispersion curves of the guided waves, obtained by Eq. ( 7) and using our algorithm, are first validated for the reference multilayer by comparing them with simulation results obtained on the Comsol ® finite element (FE) calculation software.

Figure 5 shows that a perfect agreement is obtained between the two calculation methods. This proves the reliability of the representation of the dispersion curves of the guided modes (in the case of absence of porosity, i.e. -= 0) used as a reference for all comparisons that allow to study the influence of the properties of the interface layer on the resonance modes of the structure of interest in our work.

Figure 6 illustrates the superposition of the resonance modes obtained for the reference multilayer immersed in water and the dispersion curves of the free guided waves calculated for the same structure. The resonance modes derived from the reflection coefficients, calculated with the transfer matrix method, are represented in this paper for a maximum incidence angle of 45° and a maximum frequency of 5 MHz. This incidence angle has been chosen because above this value, starting from 0.5 MHz, the reflection coefficient is equal to one and therefore no resonance modes are generated. Concerning the frequency limit of 5 MHz, it has been chosen in order to represent a sufficient number of resonance modes, with respect to the thickness of the studied multilayered plate. The global transfer matrix method suffers from numerical instabilities when the frequency-thickness product gets too high. For the structure described in Figure 4, the maximum frequency-thickness product possible before the appearance of numerical instabilities is around 39.1 MHz.mm. For a frequency of 5 MHz, the frequency-thickness product, related to the studied structure, is equal to 8.5 MHz.mm and is therefore way under the limit generating numerical instabilities.

From now on, the expression "dispersion curves of the guided waves" will refer systematically to the guided waves calculated for the reference multilayer considered in vacuum. These dispersion curves will be systematically represented by blue circles as in Figure 6. The guided modes represented in this figure are pseudo symmetric and pseudo anti-symmetric due to the absence of a plane of symmetry [START_REF] Lee | Measuring Lamb wave dispersion curves of a bi-layered plate and its application on material characterization of coating[END_REF] along the x axis. In order to identify each guided mode propagating in the multilayered reference structure, these different modes will be identified by the symbol M ‰ with 1 < i < 16. The index i increases with the frequency of appearance of the guided modes. The black patterns in this figure correspond to the minima of the reflectance function module and the white color, on the other hand, corresponds to a reflectance function modulus equal to 1. It can be seen that the minima of the reflectance function perfectly overlap the guided modes. This is in accordance with the results found in the literature [START_REF] Lee | Measuring Lamb wave dispersion curves of a bi-layered plate and its application on material characterization of coating[END_REF][START_REF] Rokhlin | Recent advances in waves in layered media[END_REF][START_REF] Lowe | The applicability of plate wave techniques for the inspection of adhesive and diffusion bonded joints[END_REF], where, in most cases, guided waves correspond to the minima of the reflectance function when the acoustic impedance of the coupling fluid is significantly lower than that of the layers forming the solid structure.

It can also be noted in Figure 6 that some guided modes, such as M and M 3 , do not have their equivalents in resonance modes. Indeed, these guided waves corresponding to horizontally polarized transverse waves [START_REF] Maghsoodi | Calculation of wave dispersion curves in multilayered composite-metal plates[END_REF], called SH (Shear Horizontal) waves, do not appear on the resonance modes since the last ones correspond to the poles of the reflectance function (Eq. 8) which was obtained considering an excitation wave polarized in the plane (x , x ), as shown in Figure 1.

Contribution of the different layers to the dispersion curves

We compare the dispersion curves of the guided waves in the reference multilayer on the one hand with those propagating in the nickel layer only (Figure 7) and, on the other hand, with those propagating in the copper one only (Figure 8). The purpose of this comparison is to identify the guided modes corresponding to waves propagating, over certain frequency ranges, only in the nickel or copper layer, or resulting from a coupling between them. This procedure allows identifying the modes in certain frequency ranges most likely to be affected by changes in adhesive layer properties, thus corresponding to a coupling between the two layers. 8 illustrate well the fact that guided modes in the multilayer get close, for certain frequency ranges, to Lamb modes propagating in the nickel and copper layers and that other modes do not correspond to any particular layer and therefore result from a coupling between the two metal layers via the interface one. Table 2 details the classification of the guided modes in the reference multilayer to the three categories described above. Acoustic coupling increases at low frequencies and the glued structure vibrates more like a "double thickness" plate [START_REF] Nagy | Adhesive joint characterization by leaky guided interface waves[END_REF], which explains why most of the coupling modes listed in Table 2 appear in the low frequency range. Indeed, in this frequency range, the resonance modes of the structure are more sensitive to the elastic properties of the adhesive layer [START_REF] Nagy | Adhesive joint characterization by leaky guided interface waves[END_REF]. In Table 2, all guided modes are analyzed with the aim of extracting the frequency ranges for which they can be considered as belonging to one of the three categories previously defined. However, guided modes corresponding to SH waves (M , M 3 , M Š , M , M , M ' ) cannot be used for adhesion characterization because the latter cannot be excited in the context of the V(z) method.

Multilayer modes

Nickel layer modes and their frequency ranges (MHz)

Copper layer modes and their frequency ranges (MHz)

Coupling modes and their frequency ranges (MHz)

‹ ! " ! ∶ $. • < • < !. ! !. ! < • < !. •' ‹ S " S ∶ $ < • < $. • ' S ∶ ' < • < z $. • < • < ' ‹ T " T ∶ $ < • < $. • " ! ∶ !. S < • < !. •z " T ∶ !. " < • < S. z ' T ∶ !. • < • < !. " C n ∶ 4.6 < -< 4.8 $. • < • < !. S S. z < • < '. • '. • < • < z ‹ ' " S ∶ S. ' < • < z ! < • < S. ' ‹ z " T ∶ !. z < • < !. "S " ! ∶ !. "z < • < S " T ∶ S. • < • < T. ' " ! ∶ T. zz < • < z ' S ∶ !. " < • < S. ' ! < • < !. z 2.4 < -< 2.8 ‹ • " z ∶ S. " < • < '. " !. 'S < • < S. " ‹ " " • ∶ S. S < • < S. z " ' ∶ T. " < • < z ' • ∶ '. "" < • < z !. " < • < S. S S. z < • < T. " ‹ • !. "• < • < z ‹ " " • ∶ T < • < z S. z' < • < T ‹ !$ " " ∶ T. $" < • < z ‹ !! " • ∶ T. TS < • < T. • " " ∶ T. " < • < z ‹ !S " • ∶ T. " < • < '. ! " z ∶ '. " < • < z '. ! < • < '. "

Influence of different parameters of the interface layer on the resonance modes

In this section, we will study the influence of different parameters, characterizing the interface layer, on the resonance modes of the multilayer. The influence of these parameters will be studied through comparisons between the corresponding modified resonance modes and the dispersion curves of the guided modes obtained in the case of the reference structure with a perfect adhesive layer. Particular attention will be paid to the propagation modes corresponding to coupling modes that are more sensitive to interface conditions. When h varies, we observe that the resonance modes are shifted over the major parts of the frequency ranges over which they are excited. For a lower value of h, the resonance modes are shifted to the higher frequencies and are generated, at a given frequency, for smaller incidence angles. This behavior means that the waves guided in the structure propagate at higher speeds as the interface layer gets thinner. Indeed, as it has been showed through Eq. ( 9), a smaller incidence angle means a higher phase velocity. The opposite behavior is observed when h increases. This might be due to the fact that, as the elastic properties of the adhesive layer tend to decrease the guided waves velocities of the structure, a thinner adhesive layer has less impact on the phase velocities of the bilayer.

Influence of the thickness of the interface layer

‹ !T " • ∶ '. ST < • < z T. • < • < '. ST ‹ !' '. '" < • < z ‹ !z '. " < • < z ‹ !• '. " < • < z
It can be seen that the variation in the thickness of the interface layer leads to a greater shift in the resonance modes corresponding to global coupling modes of the structure, such as M n (1 MHz < f < 1.5 MHz), M ' (1.42 MHz < f < 2.6 MHz), M š (1.96 MHz < f < 4.6 MHz) and M 3 (4.5 MHz < f < 5 MHz). Thus, these modes will be of great interest to be detected, in order to quantify the thickness variations of the interface layer.

As the thickness of the interfacial layer increases, some resonance modes are attenuated as illustrated by the modes M š (1.9 MHz < f < 2.56), M › (2.54 MHz < f < 3 MHz) and M (3.25 MHz < f < 4.25 MHz). This attenuation is visible through the black patterns representing the resonance modes, that transforms to grey up to light grey for some frequency ranges. This lighter color represents the fact that the reflection coefficient is higher, which is equivalent to say that more of the energy needed to excite the corresponding resonance mode is leaked inside the incident medium, thus leading to an attenuation of the mode. This attenuation may be due to the increased contribution of the interface layer which results in a higher reflectance function due to the large difference in acoustic impedance between the interfacial material and the metals constituting the other layers.

Finally, the shift in resonance modes increases with frequency, as illustrated by the M 3 coupling mode (4.49 MHz < f < 5 MHz). This is most likely due to the fact that the decrease of the wavelengths when frequency increases makes the guided modes in the multilayer become more sensitive to the interface layer.

Influence of the mass density of the interface layer

The aim here is to determine the influence of the density of the interfacial material on the propagation of guided waves in the multilayered structure. Figure 10 shows a comparison between the dispersion curves of the guided waves in the reference multilayered structure and the resonance modes obtained in the case of multilayers containing interface layers of different mass densities. It can be noted in this figure that the increase or decrease in mass density does not seem to affect the resonance modes of the multilayer. This may be due to the large difference in mass density between the nickel (ρ = 8900 kg/m ) and copper (ρ = 8936 kg/m ) layers and the epoxy (ρ = 1200 kg/m ) layer which is the matrix material of the interface layer. Thus, it appears that in the case of metallic multilayers whose effective layers all have mass densities much higher than that of the adhesive layer, it is difficult to deduce variations in the latter's density through the observation of the resonance modes of the investigated structure. 2). Since the resonance modes correspond to the minima of the reflectance function, the attenuation of the previous resonance modes, as the porosity increases, is explained by a higher amplitude of the received (reflected) signal [START_REF] Guo | The non-destructive assessment of porosity in composite repairs[END_REF]. Thus, an interface layer with greater porosity has an acoustic impedance further away from those of the metallic layers and therefore a smaller fraction of the ultrasonic energy is transmitted through the interface [START_REF] Rokhlin | Ultrasonic wave interaction with a thin anisotropic layer between two anisotropic solids: Exact and asymptotic-boundary-condition methods[END_REF]. Moreover, as porosity increases, all resonance modes are shifted to low frequencies, except for the resonance mode corresponding to M (0.6 MHz < f < 0.8 MHz), especially in the case of modes generated for higher frequencies, probably due to a higher sensitivity of the waves to the porosity of the interface layer. This behavior can be used effectively to monitor the state of an adhesive layer because, as Wang [START_REF] Wang | Evaluation of interfacial properties in adhesive joints of aluminum alloys using angle-beam ultrasonic spectroscopy[END_REF] has shown, the frequencies of occurrence of the reflectance function minima vary with the condition of the interface. Indeed, a degradation of the interface is accompanied by a strong shift of the minima in the decreasing frequency direction [START_REF] Rokhlin | Ultrasonic evaluation of environmental durability of adhesive joints[END_REF]. In addition to this shift of minima towards low frequencies, we see that resonance modes tend to appear at higher incidence angles, for a given frequency, as porosity increases, which means that their phase velocities are lower. The modes most affected by this decrease in phase velocity are M ' (1.42 MHz < f < 2.46 MHz) and M 3 (4.49 MHz < f < 5 MHz) coupling modes, sensitive to interface layer properties. Since the increase in porosity reduces the mass density of the interface layer and its shear moduli (Figure 3), it can be assumed that the decrease in resonance mode velocities is related to the decrease in propagation velocities in the interfacial material.

Thus, shift of resonance modes towards low frequencies, attenuation of coupling modes and decrease of their phase velocities can be used as indicators of the presence of porosity within the layer, corresponding to poorer adhesion.

Influence of the orientation of the pores

In the previous section, the influence of the porosity factor of the interface layer on resonance modes was presented. Now, we will vary the azimuthal angle φ that the cylindrical pores form with the incidence plane and analyze the consequences of this variation on the resonance modes. Figure 12 shows a comparison between the dispersion curves of the guided waves and the resonance modes calculated with interface layers having a porosity factor equal to 0.5 and different angles φ. A porosity of 0.5 was chosen here because this value offers a good compromise between a significant influence of the porosity of the interface layer and a reasonable attenuation of the resonance modes allowing to clearly distinguish them. Even though this porosity value might not be very frequent in real samples, except for high delamination cases, we chose it here as a case study in order to better analyze the influence of the orientation of the pores on the resonance modes. Since the interfacial material is orthotropic, it has two planes of symmetry: one with respect to φ = 0° and another with respect to φ = 90°. Therefore, we will only look at the angles φ between 0° and 90° in order to cover all cases. When φ is different from 0, we observe a decoupling of some resonance modes, such as M n at f = 1.23 MHz and f = 1.48 MHz, M ' at f = 1.89 MHz and M at f = 3.47 MHz. M n and M ' modes, for the frequencies mentioned, lie within frequency ranges assigned to the global coupling between the different layers forming the structure (Table 2). Areas where decoupling is clearly visible are highlighted in Figure 12-b). These decouplings gradually increase until φ = 45° (Figure 12-d)), then the modes resulting from these decouplings tend to be attenuated as visible on figures 12-b) to 12-e), and specifically illustrated on Figure 12e) representing the resonance modes of the structure for ž = 75°. The attenuation of the decoupling phenomenon itself, visible for ž > 45°, is also illustrated in Figure 12-f). Thus, we can affirm that the effect of the pore orientation on resonance modes is maximum for ž = 45°, and that it tends to decrease as the symmetry axes approach, which here correspond to the angular values ž equal to 0° and 90°.

It can be noted that these so-called "decoupling" modes correspond to quasi-SH waves, since they join up the SH waves velocities as it can be seen on Figure 12-d). These resonance modes do not appear among the resonance modes obtained for pores directed along the incidence plane. These modes correspond to quasi-SH waves and not purely SH waves because, as explained above, the last ones cannot be generated due to the nature of the incident wave. Thus, the fact that the cylindrical pores are offset from the incidence plane allows the generation of waves with polarizations that are almost perpendicular to the incidence plane.

The increase in φ also leads to an attenuation of the resonance modes visible through shallower minima. When φ = 0, the cylindrical pores are parallel to the incidence plane and it is possible that this configuration allows to guide more acoustic energy along the interface and therefore to generate deeper minima. To illustrate this phenomenon, the most affected resonance modes have been surrounded in Figures 12-a) and 12-f), relative respectively to the lowest φ angle (φ = 0°) and the highest φ angle (φ = 90°).

To summarize, the variation of the angle φ has consequences on the coupling of resonance modes, on their attenuation, as well as on their phase velocities. In a practical case, in order to detect a defect with a preferential orientation inside the adhesive layer, this defect will be modeled by a privileged φ orientation of the pores and the decoupling of certain resonance modes, observed experimentally, could constitute a tool to characterize this type of imperfection.

Conclusion

In this work, the effect of certain properties of an interface layer is studied in the context of an anisotropic multilayered structure. Reflectance functions and free guided waves dispersion curves are modeled using the global transfer matrix method. The characterization of variable properties of the adhesive layer by means of a comparison of resonance modes (given by the minima of the modulus of the reflectance function) with free guided waves dispersion curves of the multilayer containing a perfect epoxy layer is proposed. An anisotropic multilayered system consisting of nickel, epoxy resin and copper is the configuration studied here. The equivalence between the reflectance function magnitude minima and the free guided waves dispersion curves has been established for the considered structure. The coupling modes of the multilayer have been highlighted through a comparison between its free guided waves dispersion curves and the Lamb waves of the nickel plate only and copper plate only respectively. Then, we studied the influence of some parameters responsible for the alteration of an adhesive layer. Indeed, the impact of the thickness and density of the adhesive layer has been studied with respect to the guided modes, especially with regard to the coupling modes that are the most sensitive to the changes in the interface layer. A model assessing the interface quality through a fictitious porous material, whose cylindrical pores can be shifted with respect to the incidence plane, has also been implemented. The increase in porosity induces a shift of the resonance modes with respect to frequency and a change in the velocity of the guided modes. In addition, the variation of the orientation angle of the pores with respect to the incidence angle has consequences on the coupling of the resonance modes, on their attenuation as well as on their phase velocities. All of these parameters allow to characterize the interface layer and monitor its changes. Indeed, this model allows to model defects such as delamination and porosity in the adhesive layer, with the possibility to vary their orientation which enables us to simulate defects with a privileged orientation that can be detected through the decoupling of certain modes.

This study provides an overview of the impact of the adhesive interface layer nature on the reflectance function magnitude and the potential of the comparison with the free guided waves dispersion curves in order to characterize the interface quality within an anisotropic multilayered structure.

Thus, it could be envisaged to characterize interface defects through experimental V(z) measurements. Indeed, the reflectance functions, obtained by modeling, could be compared to the experimental ones rebuilt thanks to the inversion of measured V(z) data, a study which is expected to be performed in the near future. where:

¤ 2 = 3 1 4) 2 A.4

Figure 1 :

 1 Figure 1: Multilayered structure geometry

Figure 2 :

 2 Figure 2: Geometrical model of the composite material containing cylindrical pores[START_REF] Christensen | Mechanics of Composite Materials[END_REF]

  where M, O and C are detailed in Appendix A. Once the moduli ( , ) , * , + and + are obtained, the elastic constants & , & , & , & and & '' are calculated via the following equations: & = ( / 4* ) (5.a) & = 2* ) (5.b)

Figure 3 :

 3 Figure 3: Evolution of the shear moduli R !S and R ST with respect to the porosity

  n m ' m 3 m 3 m 3 m 33 m 3n m 3' m n m n m n m n3 m nn m n' m ' m ' m ' m '3 m 'n m ''

Figure 4 :

 4 Figure 4: Anisotropic reference multilayered structure: nickel (t u = !. S ww)/ interface layer (t x = S Rw )/ copper (t y = $. z ww)

Figure 5 :

 5 Figure 5: Comparison between the dispersion curves obtained by FE simulation and the algorithm for the reference multilayer

Figure 6 :

 6 Figure 6: Comparison between resonance modes and dispersion curves of the guided waves

Figure 7 :Figure 8 :

 78 Figure 7: Superposition of the dispersion curves of the guided waves propagating in the reference multilayer and the nickel layer

Figure 7

 7 Figure7and Figure8illustrate well the fact that guided modes in the multilayer get close, for certain frequency ranges, to Lamb modes propagating in the nickel and copper layers and that other modes do not correspond to any particular layer and therefore result from a coupling between the two metal layers via the interface one. Table2details the classification of the guided modes in the reference multilayer to the three categories described above. Acoustic coupling increases at low frequencies and the

Figure 9

 9 Figure 9 illustrates a comparison between the dispersion curves of the guided waves of the reference structure and the resonance modes obtained for different values of the interface layer thickness. The chosen thicknesses are: ™ = ™/2 and ™ = 2™.

Figure 9 :

 9 Figure 9: Comparison between the dispersion curves of the guided waves and the resonance modes calculated with different thicknesses of the interface layer: a) t = t/S and b) t′′ = St

Figure 10 :

 10 Figure 10: Comparison between the dispersion curves of the guided waves and the resonance modes calculated with different densities of the interface layer: a) oe = oe/S and b) oe′′ = Soe

Figure 11 :

 11 Figure 11: Comparison between the dispersion curves of the guided waves and the resonance modes calculated with different porosity factors of the interface layer: a) P = 0.01 ; b) P = 0.1 ; c) P = 0.3 ; d) P = 0.6

Figure 12 :

 12 Figure 12: Comparison between the dispersion curves of the guided waves and the resonance modes calculated with interface layers modeled with P=0.5 and different values of •: a) • = 0°; b) • = 10°; c) • = 30°; d) • = 45°; e) • = 75°; f) • = 90°

Table 1 : mechanical properties of the materials constituting the multilayer Material

 1

	Nickel [11]	8900	270	170	123
	Copper [20]	8936	168.4	121.4	75.2
	Epoxy [19]	1200	5.2	2.8	1.2

Mass density (kg/m ) c (GPa) c (GPa) c 33 (GPa)

Table 2 : Identification of the different categories of guided waves propagating in the multilayered structure

 2

Acknowledgement

This work has been supported by the co-labeled EURIPIDES² and CATRENE project SAM 3 (Smart Analysis Methods for 3D Integration).


Prof Laurent Fedi
email: laurent.fedi@kedgebs.com

Dr Olivier Faury
email: o.faury@em-normandie.fr

Dr Laurent Etienne
email: laurent.etienne@isen-ouest.yncrea.fr

Mapping and analysis of maritime accidents in the Russian Arctic through the lens of the Polar Code and POLARIS System

Keywords: Arctic, risk analysis, risk mapping, maritime claims, Polar Code, POLARIS

Over the last decade, more than 500 maritime accidents have been reported in the Arctic, and the political will of the Russian government to develop the Northern Sea Route (NSR) for the exportation of raw materials may lead to a rise of such events. However, data on insurance claims in the Russian Arctic are rather sketchy and not up-dated which has hindered accurate analysis so far.

Based on empirical data collection throughout the period 2004-2017, the study gathers, classifies and maps certain 'marine casualties' and 'incidents' that occurred along the NSR until the Polar Code enforcement. Three maritime events are scrutinized through the lens of the Polar Code risk mitigation tools and POLARIS system in order to determine if the joint use of these new prophylactic instruments would have prevented their occurrence.

As preliminary results, the research identifies the primary causes of identified maritime accidents that occurred in the Russian Arctic and provides their categorization in accordance with IMO standards. It appears that serious 'marine casualties' represent the greatest number of insurance claims. Shedding light on the type of risk occurrence a ship may face along the NSR, it brings valuable implications for maritime operators who are still awaiting historical data for more realistic risk assessment. Finally, it demonstrates that the combination of the Polar Code and POLARIS system is an appropriate tool for risk mitigation and encourages stakeholders to implement best practices for Arctic shipping in their common interest.

Introduction

Due to global warming, sea ice coverage and ice thickness in the Arctic Ocean have progressively decreased in recent decades (IPCC, 2018;[START_REF] Serreze | Arctic sea ice trends, variability and implications for seasonal ice forecasting[END_REF]Corbett et al., 2010). This offers new economic opportunities in particular for the shipping industry [START_REF] Lasserre | Compared transit traffic analysis along the NSR and the NWP[END_REF]Eguíliz et al., 2016;[START_REF] Melia | Sea ice decline and 21 st century trans-Arctic shipping routes[END_REF][START_REF] Kum | A root cause analysis for Arctic marine accidents from 1993 to 2011[END_REF] both for the Northwest Passage (NWP) and the Northern Sea Route (NSR). The construction of a new Liquid Natural Gas (LNG) plant in the Yamal Peninsula (Russia), allowed the first operations through the Bering Strait in July 2018. Furthermore, even if the target of the Russian Federation was to export between around 130 and 150 million tons by 2030 (Ministry of Energy of the Russian [START_REF] Federation | Energy Strategy of Russia for the period up to 2030[END_REF], figures were revised to 80 million tons by 2024 [START_REF] Staalesen | Shipping on Northern Sea Route up 40%[END_REF] which is still considered as highly ambitious (Staalesen, 2019a). The NSR gives the opportunity to reduce distances [START_REF] Zhang | Is the Northern Sea Route attractive to shipping companies? Some insights from recent ship traffic data[END_REF], and fuel consumption thereby lowering carbon emissions [START_REF] Kitagawa | Arctic routing: challenges and opportunities[END_REF]. Although the use of the NSR is questioning notably for international transit between non-Russian ports in Europe and Asia (Lasserre, 2016: Faury and[START_REF] Faury | The Northern Sea Route competitiveness for oil tankers[END_REF]Cariou and Faury, 2019), this route crystallizes both greater maritime actors' interests, among which oil and gas producers that represent the lion's share of the goods shipped along the NSR but also cruise companies [START_REF] Cajaiba-Santana | The emerging cruise shipping industry in the arctic: Institutional pressures and institutional voids[END_REF] and container shipping lines such as Maersk and Cosco.

Nevertheless, a high level of risk in particular as regards ice characterizes the NSR and the Russian Arctic as a whole (ABS, 2014) despite limited statistics on shipping accidents (AMSA, 2009;[START_REF] Goerlandt | An analysis of wintertime navigational accidents in the Northern Baltic Sea[END_REF]).Yet, the number of insurance claims doubled in Arctic Circle waters throughout the period 2007-2017[START_REF] Allianz | Safety and shipping review: an annual review of trends and developments in shipping losses and safety[END_REF] thus highlighting an opposite trend compared to global safety improvements justifying the adoption by the International Maritime Organization (IMO) of the Polar Code (PC) (IMO, 2014a) designed to reduce operational risks in Polar waters via safety and environmental prevention measures (Fedi et al., 2018;[START_REF] Dalaklis | The future of Arctic shipping business and the positive influence of the international code for ships operating in polar waters[END_REF].

The aim of this paper is twofold. First, it carries out an in-depth analysis of certain maritime accidents that occurred in the Russian Arctic, identifying their root causes, level of seriousness and locations. The analysis enables a taxonomy of these events in accordance with the Casualty Investigation Code [START_REF] Imo | Resolution MSC.255[END_REF]. It seems that existing literature does not provide this categorization. Secondly, three specific maritime claims linked to the ice risk are analyzed through the lens of the "Polar Operational Limit Assessment Risk Indexing System" (POLARIS) in order to verify if their occurrence would have been prevented. Even though IMO promotes POLARIS as well as classification societies, so far, little research has been carried out on this valuable system for decision-making of Polar navigation. The contribution of this paper is to highlight the relevance of the combination of PC provisions with POLARIS instrument in formal risk mitigation for the NSR and particularly in ice-covered areas.

After this introduction, Section 2 reviews current literature on Arctic shipping risks with a particular focus on the NSR. Section 3 deals with the methodology. The analysis of maritime accidents occurring in the Russian Arctic and the case studies are contemplated in Section 4 while Section 5 provides a discussion and managerial implications. Finally, some concluding remarks are addressed in Section 6.

Literature review

As illustrated by the recent adoption of the PC, the safe navigation in Arctic waters has become a pressing issue of concern. The ice presence, remoteness [START_REF] Haavik | Remoteness and Sensework in Harsh Environment[END_REF]Lasserre et al., 2016), lack of search and rescue (SAR) infrastructures (SARC, 2017), low temperatures, drifting icebergs -over 1,000 icebergs drifted into North Atlantic shipping lanes in 2017 (ALLIANZ, 2018) -, rapidly changing and severe weather conditions and extended period of darkness represent considerable hazards for ships complicated by the lack of experience of crew [START_REF] Fu | Towards a Probabilistic Model for Predicting Ship Besetting ice in Arctic Waters[END_REF][START_REF] Montewka | Towards a Probabilistic Model of a Ship Performance in Dynamic Ice[END_REF][START_REF] Tikka | Tanker Design Considerations for Safety and Environmental Protection of Arctic Waters: Learning from Past Experience[END_REF]. Moreover, as predicted by some scholars [START_REF] Borgerson | Arctic meltdown: the economic and security implications of global warming[END_REF][START_REF] Balto | Risk of Ship Collision in the Barents Sea in 2030[END_REF], the number of maritime accidents has recently increased. Although in 2006 there were only eight, they reached 71 in 2017 [START_REF] Allianz | Safety and shipping review: an annual review of trends and developments in shipping losses and safety[END_REF], meaning that accidents skyrocketed within one-decade contrary to the long-term downward trend of worldwide losses [START_REF] Iumi | Global marine insurance report[END_REF][START_REF] Allianz | Safety and shipping review: an annual review of trends and developments in shipping losses and safety[END_REF]. Accordingly, this specifically raises the question for the NSR in light of the hazardous goods carried (oil and gas) and a greater traffic likelihood in coming years.

Russian Arctic shipping risks and accidentology

Located in the Russian Federation Exclusive Economic Zone (EEZ), the NSR represents several navigational routes. Pursuant to the Northern Sea Route Administration (NSRA, 2012), the NSR is "a water area adjoining the northern coast of the Russian Federation, including internal sea waters, territorial sea, contiguous zone and EEZ […], and limited in the East by the line delimitating the sea areas with the USA and by the parallel of the Dezhnev Cape in the Bering Strait; in the West, by the meridian of the Cape Zhelanie to the Novaya Zemlya archipelago, by the east coastal line of the Novaya Zemlya archipelago and the western limits of the Matochkin Shar, Kara Gates, Yugorski Shar Straits". The NSR is described as running through four main zones: the Kara, Laptev, East-Siberian and Chukchi Seas. The NSR extends for around 3,000 miles and its navigation must comply with Russian legislation. Connecting Asian and European ports, annual navigation along the NSR is feasible depending on the ice or polar class of the vessel. From 2010 to 2017, 302 vessels performed along the NSR between European and Asian ports, between Russian and European or Asian ports and between Arctic ports themselves with a significant 40% increase in 2017 (NSRA, 2018). The main ports in the Russian Arctic saw a 40% rise in goods handled in 2016 compared to 2015 [START_REF] Sarc | Arctic Search and Rescue Capabilities Survey[END_REF]. In 2018, roughly 20 million tons [START_REF] Staalesen | Shipping on Northern Sea Route up 63 percent[END_REF] passed through Russian Arctic seaports and was made up mostly of LNG. Thanks to the opening of production at the Yamal LNG plant, traffic is forecasted to grow even more rapidly over the coming years [START_REF] Iudin | The Russian Sea Areas and Activity Level Up to[END_REF] while container shipping lines show little interest in this route (Lasserre, 2016).

The harsh and variable conditions along the NSR "challenge mariners, technology and systems" even when the NSR is open (AMSA, 2009;[START_REF] Abs | Navigating the Northern Sea Route -Status and Guidance[END_REF]. The main factor influencing the choice of the NSR is the distribution of ice [START_REF] Liu | The potential economic viability of using the Northern Sea Route as an alternative route between Asia and Europe[END_REF] and ice may be encountered at any time including summer where ice massifs are frequent [START_REF] Abs | Navigating the Northern Sea Route -Status and Guidance[END_REF]. Despite a longer seafaring season, waters are infested by floating ice and ice conditions are still extremely challenging for crew members, vessels (Marchenko, 2012a) and make the seafaring conditions somewhat unpredictable [START_REF] Abbassi | Risk analysis of offshore transportation accident in Arctic waters[END_REF]. The ice floes or ice drift may put a considerable pressure against the ships' hulls. Vessels must adapt their route, trying to sail a way of "least resistance through the ice mass" to avoid the risk of being stuck or limit any risk of incident whatsoever, let alone accident [START_REF] Vihanninjoki | Arctic Shipping Emission in the Changing Climate[END_REF]. Consequently, officers face difficulties to stick to the planned route chosen before the vessel's departure. Additionally, the lack of infrastructures and the remoteness of the NSR render environmental conditions more challenging for seafarers, especially in case of unexpected events (Marchenko, 2014a;[START_REF][END_REF]. The NSR counts only eight major ports [START_REF] Gritsenko | Policy environment analysis for Arctic seaport development: the case of Sabetta (Russia)[END_REF] subject to hard conditions and entering or operating in their facilities is often impossible.

On another note, few weighty studies reported the root causes of marine accidents in Arctic waters and in particular in the Russian Arctic. Although the Arctic Marine Shipping Assessment (AMSA) report identified more than 290 accidents over the period 1995-2004 (2009), there was a 'complete absence' of reported events in the Russian Arctic. The AMSA thereby recommended the development of an accurate database. Notwithstanding the existence of the IMO GISIS, the European EMCIP or the Lloyd's Marine Intelligence Unit database, numerous 'weaknesses' characterize these databases [START_REF] Luo | Half-century research developments in maritime accidents: Future directions[END_REF] on the Arctic region due to underreporting [START_REF] Goerlandt | An analysis of wintertime navigational accidents in the Northern Baltic Sea[END_REF] as well as different information and taxonomy applied to the same maritime events [START_REF] Ladan | Data Sources for Quantitative Marine Traffic Accident Modeling[END_REF]. According to [START_REF] Goerlandt | An analysis of wintertime navigational accidents in the Northern Baltic Sea[END_REF] limited and incorrect data lead to gaps in the analysis and may distort results. Furthermore, these databases are not often publicly accessible. Hence, providing an accurate statement on an accident in a given area is still a difficult task [START_REF] Grabowski | Human and organizational error data challenges in complex, large-scale systems[END_REF].

With the exception of [START_REF] Marchenko | Russian Arctic Seas: navigational conditions and accidents[END_REF]2015) who listed 90 accidents in the Russian Arctic between 1990 and 2010, few studies have been carried out up to now in particular with a global approach including identification, categorization and mapping of claims that occurred along the NSR. Most current reports use the notions of 'incident', 'accident' and 'casualty' alternatively and concurrently. This lack of data combined with approximate classification constitutes a significant gap in the existing literature insofar as it hinders stakeholders from precisely assessing the actual level of accidentology in the Russian Arctic and then implementing adapted responses. It is recognized that this sketchy statistical data complicates the shipowners' and underwriters' decision-making (Fedi et al., 2018a;[START_REF] Kiiski | Feasibility of Commercial Cargo Shipping along the Northern Sea Route[END_REF][START_REF] Sarrabezoles | Arctic shipping insurance: towards a harmonisation of practices and costs?[END_REF] on appropriate price quotations and operations.

The prophylactic responses: the PC and the POLARIS system

The new PC aims to respond to increased maritime activities in Polar Regions (Fedi et al., 2018;[START_REF] Dalaklis | The future of Arctic shipping business and the positive influence of the international code for ships operating in polar waters[END_REF][START_REF] Fedi | Les Principaux Enjeux et Impacts du Code Polaire OMI[END_REF]. Following a 'risk-based approach', the code pinpoints the main hazards in these zones, their likely impacts and provides appropriate measures to mitigate them in a prophylactic way (Fedi et al., 2018a;2018b;[START_REF] Fedi | Arctic Shipping Law From Atomised Legislations to Integrated Regulatory Framework: The Polar Code (R)Evolution?[END_REF]. Fedi, Faury and Gritsenko (2018a) demonstrated that the PC eased risk appraisal and could be a response to lack of data through mandatory safety and environmental standards. Defined as a 'toolbox' reinforcing polarseaworthiness (Cullen, 2015), they demonstrated that the PC, in codifying risk assessment, represented a key determinant of marine insurance through 'proceduralization of risk'. This proceduralization lies in three main tools (PC, Part I-A): the 'Polar Ship Certificate' (PSC) attesting the ship's ice class (A, B, C) depending on her capabilities to operate in polar waters and the 'Polar Water Operational Manual' (PWOM) establishing detailed safety requirements pursuant to design capabilities and limitations that the master and crew must comply with in order to take "sound operational decisions" (DNV-GL, 2017). Finally, 'voyage planning' must be defined that encompasses nine requirement categories including hydrographic information, extent and type of sea ice, presence of icebergs, temperatures, aids to navigation, places of refuge or SAR capabilities along the chosen route.

In conjunction with the PC, one observes the POLARIS system developed by the International Association of Classification Societies (IACS) and some Arctic countries (IMO, 2014b(IMO, , 2014c) that defines a methodology aiming at evaluating the ship's operational capabilities in ice [START_REF] Imo | Guidance on Methodologies for Assessing Operational Capabilities and Limitations in Ice[END_REF]. POLARIS compares the different ice types to the vessel class with the aim of choosing the safest navigation route and 'the optimal class of the vessel willing to sail within the polar waters' (Fedi et al., 2018b). Based on ice conditions, POLARIS may be useful for classification and insurance companies according to [START_REF] Kujala | Challenges for Application of Risk Based Design Approaches for Arctic and Antarctic Operations[END_REF] and Fedi et al. (2018b). For the latter, POLARIS constitutes "a fairly complex and multipurpose tool" both upstream and downstream of the shipowner's decision process for safer Arctic navigation. Building on this basis, it is assumed that POLARIS system combined with the PC's tools can reduce risk occurrence along the NSR where ice is predominant and would have prevented some accidents from occurring prior to the complete enforcement of the PC.

Methodology

This paper takes a multi-method approach with a qualitative and quantitative research design. The qualitative part firstly deals with the data collection of shipping accidents in different sources such as the International Union of Marine Insurance (IUMI), accident archives especially from the Protection of the Arctic Marine Environment (PAME), marine investigation reports from CLARKSONS, IHS Markit database (2019), the IMO GISIS, a recent study on Russian Arctic claims [START_REF] Fedi | Mapping and analysis of maritime claims in the Russian Arctic based on POLARIS System[END_REF], and the trade press between 2004 and 2017. One justifies this period since the 2009 AMSA report provided a global analysis up to 2004 and the entry into force of the latest provisions of the PC was effective in July 2018 for manning and training requirements whereas safety and environmental obligations were implemented in January 2017 (DNV-GL, 2017).

We secondly categorized the accidents in accordance with the appropriate terminology based on the IMO Casualty Investigation Code (IMO, 2008) and effective from 1 st January 2010 as explained in Table 1.

Insert Table 1 here

The third stage was to map maritime events in order to provide their location along the NSR, to verify if certain areas and periods are riskier than others and what key or aggravating factors could be identified.

The quantitative part of the methodology is focused on the analysis of vessels that faced an accident or incident related to navigation in dangerous areas for which we contemplate their Risk Index Outcome (RIO) as defined by POLARIS. The RIO is a value designed to assess operating limitations in ice depending on three levels: 'normal operation', 'elevated operational risk' and 'operation subject to special consideration'. This RIO can be computed using a risk index table as defined in [START_REF] Imo | Guidance on Methodologies for Assessing Operational Capabilities and Limitations in Ice[END_REF]. When the ice thickness is above the ship ice class capability, the risk index is negative. Respectively, if the ship has the ability to sail through an ice thickness, then the risk index is positive. This risk index is then multiplied by the sea ice concentration to compute the POLARIS RIO.

To define the RIO in this paper, ice conditions statistics (thickness and concentration) were extracted over a period of 28 years on a daily basis from Copernicus database (2019) and aggregated in a grid of 12.5 km cells. As the European Union's Earth Observation Program, Copernicus provides open real-time data especially on atmosphere, land and marine environments (EU Commission, 2015). Based on [START_REF] Faury | La Route Maritime Nord : complément ou concurrent des routes conventionnelles ?[END_REF], we defined the level of risk of each case, the day of the accident and seven days before and after to determine whether the situation could have been anticipated or not. Ice conditions were evaluated in the area of the event according to three scenarios: worst case, median and best case based on POLARIS system. These three scenarios are aligned with the observed descriptive statistics (worst case = 25% lower values, median case = 50%, best case = 75% upper values) of the aggregated POLARIS RIO values within each grid cells for each day of year (1 to 365) over 28 years.

Pursuant to the IMO Guidance, three out of the five elements characterizing POLARIS deal with ice: the ice type definitions pursuant to the World Meteorological Organization (WMO), the different ice regimes such as partial ice concentrations or ice-free waters and ice decay in warmer temperatures (Fedi et al., 2018b). We considered POLARIS as a decision tool through risk classification and proceduralization (Fedi et al., 2018b) since it provides the RIO index that represents the cornerstone of the decision to fully operate or to limit shipping operations as shown in Table 2.

Insert Table 2 here

Vessels under PC1-PC7 categories (corresponding to the A and B categories of the PC) are designed for operations in polar waters under certain ice conditions (DNG-GL, 2017). For RIOs ≥ 0, there are no specific requirements, sailing is allowed as 'normal operation'. For other RIOs between ≥ 0 and -10, sailing is 'subject to special consideration' that is to say with icebreaker assistance for instance and will depend on the vessel ice class since operational risk is elevated. Finally, if the RIO is below -10, no navigation must be undertaken. However, the final decision will belong to the master and whatever the ice class of the ship (Fedi et al., 2018b).

Analysis of maritime claims occurring in the Russian Arctic

Risk categories

Sailing within the Russian Arctic implies dealing with numerous risks mostly identified within the PC and that can be organized in four main categories: climate, technical, human and geographic. Yet looking at parameters that directly impact navigation, climate appears as one of the major risks, ice having a strong impact on navigation within polar waters (Cariou et al., 2020;[START_REF] Faury | The Northern Sea Route competitiveness for oil tankers[END_REF][START_REF] Sarrabezoles | Arctic shipping insurance: towards a harmonisation of practices and costs?[END_REF]. Floating ice has been identified as the key factor of most maritime casualties [START_REF] Marchenko | Russian Arctic Seas: navigational conditions and accidents[END_REF](Marchenko, , 2014b)). Depending on ice thickness and concentration a vessel may not be able to sail [START_REF] Löptien | Ice and AIS : Ship speed data and sea ice forecast in the Baltic Sea[END_REF]. However, numerous measures exist to prevent these risks such as the use of polar-class vessels with reinforced hull and a more powerful main engine or the assistance of icebreaker made sometimes compulsory by underwriters (Fedi et al., 2018a) and NSR requirements (NSRA, 2012).

With regard to the current literature ice appears as the main risk weighing upon ships. Besides, the fast-changing conditions due to freezing and thawing of the ice blanket can position the vessel in unexpected situations. The ice extent changes on a yearly basis and the months with the smallest and largest areas of ice are respectively September and March (Figure 1). A considerable variability in the ice expanses can be noticed between May and August and between October through December while the ice extent amplitude is less extensive during the rest of the year.

Insert Figure 1 here Fedi et al. (2018a) highlighted that if the melting ice provides a larger period of navigation, it also generates an increase of ice drifting which can be a source of accidents [START_REF] Faury | Risk management in the Arctic from an underwriter's perspective[END_REF][START_REF] Dalaklis | The future of Arctic shipping business and the positive influence of the international code for ships operating in polar waters[END_REF] and thereby an unsuitable speed with such ice conditions (Marchenko, 2014b). Consequently, the ice and fast-changing climate conditions can block a vessel and impact the integrity of the hull or other parts such as propellers. Furthermore, the risks related to the dangerous and polluting nature of hydrocarbons and mineral products mainly carried along the NSR cannot be occulted in light of their potential environmental impacts in case of accident.

Data collection and mapping on maritime claims along the NSR

Based on multi-source databases already mentioned and [START_REF] Fedi | Mapping and analysis of maritime claims in the Russian Arctic based on POLARIS System[END_REF], we identified 36 events involving SOLAS vessels that occurred along the NSR from 2004 to 2017. The ship's name, its type, its ice class, year of event, and factor, that is to say the root cause of these events were collected. We classified them either as 'Marine Incident' (MI), 'Marine Casualty' (MC), 'Serious Casualty' (SC) or 'Very Serious Casualty' (VSC), depending on their respective consequences pursuant to the Casualty Investigation Code (IMO, 2018). SC represented 75% of the claims leading to 14 collisions while MI only concerned 8 events (19%) mainly involving ships stuck in ice (85%) and 2 VSC (sinking and death) were reported.

Looking at the factors, climate conditions were dominant and at the origin of 85% of MI, 44.4% of SC and 100% of VSC respectively. In most cases, ice represented 65% of claims. Ice impacts were not homogeneous leading either to MI without detrimental consequences for ship and crew such as drift, ship stuck or to SC involving grounding and collision. The other key factor was human for 22 % of the cases and responsible for SC such as injury (2), fire (2) and collision (4). Mechanical issues accounted for 13% of MI (1) and SC (4) contributing to vessel's drift (MI), fire, mechanical damage and miscellaneous issues (SC).

Insert Figure 2 here

As illustrated in Figure 2, the claims are progressively fewer in the northern part of the NSR during the period 2004-2017. The Barents Sea concentrated the highest number of claims with a total of 13 out of 36 representing 36% of events and can be considered as the most significant accidental area along the NSR. Five claims were due to climate conditions and 5 to human error. The Kara Sea can be analyzed as the second most dangerous area with 6 claims accounting for 20% of the total. This sea is made up of several straights through isles and archipelagoes creating difficult sailing conditions especially when ice is present. Climate remained the determining factor in the Kara Sea (for 3 claims out of 7). Ice was also involved in three claims in the Pechora Sea (1) and the White Sea (2). As regards the Chukchi Sea, only three claims are recorded and worth 8.3% of the total. Nevertheless, SC systematically occurred in this area, mainly collisions, ice factor contributing to two thirds of events.

The fact that most of the casualties occurred in the Western part of the NSR (24 out of 36) makes sense since the main oil and gas fields are exploited in this area. Indeed, we could cite the shuttle tankers between Prirazlomnoye in Pechora Sea and Murmansk or sailing of tankers between Novy Port (Yamal Peninsula) and Murmansk. Moreover, if navigation in convoy with an icebreaker mitigates the risk represented by ice, it may be at the origin of collisions [START_REF] Zhang | Safety distance modeling for ship escort operations in Arctic ice-covered waters[END_REF] 1 .

Case studies: evaluation analysis and discussion

As explained previously, tankers, general cargo and dry bulk represent the largest share of accidents faced by SOLAS vessels in the Arctic. We analyzed three cases: Sinegorsk (bulk carrier), Inger (general cargo), and Chukotka + (tanker). Insofar as ice and human error are the two main factors of risk occurrence, these cases are analyzed through the scope of 1 http://arctic-lio.com/wp-content/uploads/2018/05/news_2017-04-25_Helsinki.pdf POLARIS in order to determine if they could have been avoided and to evaluate the potential limits of this system as a decision tool for deck officers in light of PC provisions.

Table 3 summarizes the RIO of each vessel pursuant to its ice class and the ice conditions encountered when the accident occurred. To define the predictability of the event, we computed the RIO one week before and after the event, and the probability of ice conditions. We assessed the RIO according to worst, median and best-case scenario.

Insert Table 3 here

Case analysis 4.3.1.1 Inger case

On 19 March 2006, the general cargo Inger was sailing in the White Sea while it sustained damage following contact with ice. According to the POLARIS RIO, the Inger should have experienced favorable operating conditions as shown in Figure 3. However, its 11.83 RIO (Table 4) means that ice drifting may have existed. Indeed, the vessel collided with drifting ice that seriously impacted the hull and the propeller shaft. The ship was stopped for two months in a shipyard and returned to service on 17 May with detrimental economic consequences.

As previously explained, this event puts emphasis on the importance of voyage planning as imposed by the new PC and indirectly raises the question of officers' experience in polar navigation that is now framed by the said PC.

Insert Figure 3 here

Sinegorsk case

The M/V Sinegorsk was sailing from Arkhangelsk to Pevek in convoy with the M/V Johann Mahmastal assisted by the diesel icebreaker Admiral Makarov in order to supply cargo. While the goods were supposed to be delivered on 30 December 2016, they arrived on 7 January 2017.

The Sinegorsk is an ice-class 1AS vessel meaning that according to [START_REF] Konygin | Oil Tanker Transportation[END_REF], it shall not be able to sail in ice with a thickness above 0.8 meters. Table 3 stresses that the situation faced by the vessel could have been anticipated. In fact, one week before and after the RIO is below -10, which is its maximum limit. Moreover, based on ice conditions over the last 28 years, one notices that the area is barely navigable for an ice class 1AS in two scenarios out of three (Table 3).

Figure 4a highlights that the vessel is surrounded by ice conditions representing a direct danger. In other words the 7.44 RIO has to be put in perspective with the RIO in the vicinity. This approach is confirmed by Figure 4b and 4c that show the existence of red zone remaining close to the bay entrance. Indeed, ice conditions changed and the ice thickness reached 1m with ridges at 2.5 meters, blocking the convoy at the entrance to Chaunskaya Bay, 24 miles from Pevek2 and leading to an eight-day delivery delay. In such circumstances, if the PC had been in force, the mandatory PWOM would have required anticipated scenarios of operating and environmental conditions.

Insert Figure 4a here Insert Figure 4b,4c, 4d here

Chukotka +

On 30 October 2017, as it was sailing, the tanker Chukotka + became stuck in the ice close to the Sannikov Strait. Regarding POLARIS RIO (Figure 5), the vessel was sailing in a safe area and the ice conditions that day were median (10.10). A 45 cm ice thickness covered 99% of the area. As explained by Faury and Givry, (2017), a 1A should not be able to sail on its own with a thickness above 30 cm. Even though the tanker was not damaged nor did an oil spill occur, the vessel had to wait 12 days for the assistance of the nuclear icebreaker Yamal3 and faced related unexpected expenses.

Consequently, this claim stresses the importance of the captain's ability to assess the ship's capabilities. If the PC had been in force, the master would have defined voyage planning beforehand and respected the PWOM with related operational limitations.

Insert Figure 5

Discussion and managerial implications

Whereas the recorded shipping accidents do not pretend to constitute an exhaustive database, their analysis provides a basis on NSR accidentology pursuant to the IMO taxonomy during the 2004-2017 period. Obviously, most of these events are 'casualties' illustrating their seriousness. With regard to accident factors, our assessment highlights that climatic conditions and ice remain key determining causes of risk occurrence (AMSA, 2009;Marchenko, 2012a;[START_REF] Abs | Navigating the Northern Sea Route -Status and Guidance[END_REF]Abassi et al., 2017). Collisions represent the most frequent SC while ships stuck and drifts constitute the common MI showing that ice conditions along NSR are still challenging for ships [START_REF] Vihanninjoki | Arctic Shipping Emission in the Changing Climate[END_REF]Marchenko, 2014aMarchenko, , 2014b)). Moreover, human factor is identified as the second cause of claims in the Russian Arctic and emphasizes the importance of training and experience in Polar areas [START_REF] Tikka | Tanker Design Considerations for Safety and Environmental Protection of Arctic Waters: Learning from Past Experience[END_REF][START_REF] Kum | A root cause analysis for Arctic marine accidents from 1993 to 2011[END_REF] as currently required by the PC (Fedi et al., 2018a).

With regard to the typology of vessels, our analysis shows that 25% of events involve tankers, 19.4% general cargo and 16.7% bulk carriers. Insofar as raw materials are mainly carried along the NSR [START_REF] Faury | The Northern Sea Route competitiveness for oil tankers[END_REF], this statement is hardly surprising. Nevertheless, tankers face most SC in particular collisions. It raises the topical issue of oil spills and their environmental impact in fragile Arctic eco-systems (AMSA, 2009;2017). Considering the strategic role of the NSR in the exportation of Russian oil products, this traffic is expected to grow while related risks shall proportionally increase. A worst-case scenario cannot be excluded (Fedi et al., 2018a;[START_REF] Johannsdottir | Systemic risk of maritime-related oil spills viewed from an Arctic and insurance perspective[END_REF] and the low density of SAR infrastructure heightens the threat of oil spills [START_REF] Vihanninjoki | Arctic Shipping Emission in the Changing Climate[END_REF]. Finally, one observes the absence of passenger vessel claims along the NSR with contrast to the 27 events recorded by AMSA (2009).

Through the case studies, we underlined the potential contribution of POLARIS as a risk mitigation tool for ice-covered areas and the importance of officers' competencies and experience. POLARIS implementation prior to 2018 would certainly have prevented accidents occurring along the NSR where ice still involves huge risks. Nevertheless, as the cases implicitly revealed, this system shows some limitations and requires complementary measures for an appropriate decision-making. POLARIS is not self-sufficient since it only covers ice conditions which constitute a sole parameter of ship's operational limitations (Fedi et al., 2018b). Whereas the human factor is particularly important, it is not included in POLARIS. The AMSA report (2009) highlighted the determining role of human factor with 76.7% of the total number of accidents in Arctic waters. That is why POLARIS needs be apprehended and completed by PC tools. According to the PC, officers must indeed consider the vessel's capabilities to sail through or in ice-covered areas thanks to the PSC in particular and also evaluate a wide range of operating and environmental conditions beyond ice. This relates to low air temperature, high latitude, remoteness or weather conditions. In addition, appropriate voyage planning must be conducted and officers shall respect the PWOM and related procedures on operational limitations based on anticipated scenarios addressing key risks and normal and emergency conditions (DNV-GL, 2017). If the worst-case scenario was further envisaged by officers, it would allow different decisions such as the choice of different routes or the recourse to icebreaker assistance. Here, the PC requirements on specific polar crew's certification compensates for the POLARIS shortage on human element. Accordingly, this system cannot be 'disconnected' from the key PC provisions such as PSC, PWOM, voyage planning, experienced and trained crew members in treacherous waters still poorly or uncharted (hydrography and bathymetry) with frequent fogs, icy blizzards, and blocks of drifting ice [START_REF] Abs | Navigating the Northern Sea Route -Status and Guidance[END_REF]. Thus, the implementation of POLARIS combined with the PC enable operators to adopt a proper formal risk assessment for undertaking a journey via the NSR and to prevent risk occurrence in a prophylactic way.

Furthermore, the beneficial impacts of POLARIS system must be evaluated in the medium and long-run and especially in light of greater maritime activities in the Russian Arctic. As a modern methodology for assessing ship's operational capabilities, the benefits of POLARIS must be promoted by Arctic stakeholders, classification societies, shipowners and underwriters in particular. This promotion seems necessary insofar as POLARIS is not mandatory. As stated in the PC Recommendations Part-I-B (IMO, 2014), other methodologies can be used. The Arctic Shipping Best Practice Information Forum that aims to promote the PC and to support some form of good practice for Arctic shipping (AMSA 2017;Fedi et al., 2018a) provides a relevant setting for this promotion.

These preliminary results should allow Arctic states, maritime operators such as shipowners, classification societies and underwriters who need statistics in particular [START_REF] Kiiski | Feasibility of Commercial Cargo Shipping along the Northern Sea Route[END_REF]Fedi et al., 2018a) to have a better operational knowledge on past claims that occurred along the NSR (Schmied et al., 2016). Nevertheless, if universal regulation such as the PC and more realistic operational risk appraisal as POLARIS are necessary, they are not sufficient to holistically manage Arctic shipping risk. Through a general examination, one should not overlook potential severe accidents [START_REF] Marchenko | Maritime activity in the high North -the range of unwanted incidents and risk patterns[END_REF]Fedi et al., 2018a) and notably underestimate the dangers of some bulks (e.g. oil products and minerals) to the environment in case of accident. It is not necessary to wait for a second Exxon Valdez or Titanic before making changes. Public and private stakeholders must proactively work in close collaboration not only for prevent and reduce current casualty rates in the coming years but also to provide the required responses to emergency preparedness with adequate places of refuge and SAR infrastructures [START_REF] Gritsenko | Policy environment analysis for Arctic seaport development: the case of Sabetta (Russia)[END_REF].

Finally, pursuant to the PAME 2019-2021 Work Plan (PAME, 2019), the Compendium of Arctic Shipping Accidents (CASA) should be created with the purpose of providing up-dated data on the Arctic Ocean. It should constitute a valuable source of information for more exhaustive analysis on the nature of accidents, their location, their impacts and contemplating measures that might be pursued to reduce their occurrence. If Russia participated in this database, it would allow an accurate identification and mapping of shipping accidentology along the NSR with a greater transparency.

Concluding remarks

So far, few works of research have been carried out following the IMO Casualty Investigation Code and providing a global approach of accidents in the Russian Arctic. The analysis gave an overall picture of the situation and underlined that serious casualties represented the largest number of events. Our results confirmed the role of ice as determining factor in marine casualties during the 2004-2017 period. In light of growing traffic along the NSR, these preliminary results allow stakeholders to have a better knowledge on past accidents and to carry out a better formal risk assessment for undertaking a journey via the NSR. This risk evaluation must not be neglected considering that the NSRA has recently observed numerous violations of rules (NSRA, 2017).

The second contribution of the research was to highlight the relevance of the POLARIS instrument combined with the PC in the shipping risk mitigation for the NSR ice-covered areas. Through the case studies, we have shown that the couple POLARIS -PC constituted a valuable decision support tool and an appropriate mitigation response to ice risk which remains one of the greater hazards in the Russian Arctic. Yet we also emphasized the importance of the officers' skills and actual experience to assess situations properly and the capacity of their vessel to sail in risky areas. Moreover, sailing in apparent safe areas does not mean that there is no danger as the Inger and Chukotka + cases pointed out.

Regarding the limitations of the study, the provided information is not exhaustive and does not constitute a database while different sources were used in order to optimize the statistics. Besides, our analysis is mainly based on three case studies. Consequently, future research should evaluate if the same results can be duplicated in other case studies with similar contexts. In addition, considering the recent entry into force of the PC, operators are in a transition phase and a certain time is necessary to assess the long-term impacts on the NSR. Nevertheless, the new prophylactic PC provisions associated with POLARIS should contribute to a reduction of insurance claims in the long run. This analysis is a part of our research agenda and should be facilitated by the access to the promising Compendium of Arctic Shipping Accidents.

Figure 1 :

 1 Figure 1: Arctic ice extent in 1996, 2006 and 2012

Figure 2 .Figure

 2 Figure 2. Accident representation along the Northern Sea Route

Table 1 . Taxonomy of marine incidents and casualties Taxonomy Marine incident Marine casualty

 1 Authors based on IMO Resolution MSC.255(84) 16 May 2008

			Serious casualty	Very serious casualty
	Definition -Endangers the	-Fire	-Total loss of the ship
	safety of the ship,	-Explosion	-Death
	its occupants or	-Collision	-Severe damage to the
	any other person	-Grounding	environment
	or	the	-Contact	-Loss, presumed loss or
	environment		-Heavy weather damage	abandonment of a ship
			-Ice damage
			-Hull cracking
			-Suspected hull defect
			-Material damage to a
			ship
			-Stranding or disabling of
			a ship
			-Material damage to
			marine infrastructure
			-Severe damage to the
			environment
	Source:

Table 2 . Risk Index Outcome Criteria

 2

	Source: IMO (2016)

https://gcaptain.com/russian-vessel-convoy-beset-ice-northeast-passage/

https://thebarentsobserver.com/en/arctic/2017/12/emergency-arctic-waters-oil-tanker-gets-stuck-ice-driftsashore