Numerical prediction of ice rubble field loads on the Norströmsgrund lighthouse using cohesive element formulation

An attempt has been made to predict the ice rubble field load on Norströmsgrund lighthouse by using Cohesive Element (CE) formulation. Two sub-load events were selected to validate the numerical and material model used in simulation of interaction of the ice rubble field and lighthouse. A literature review of simulation of rubble field structure interaction methods is also included in order to illustrate the knowledge gaps and highlight shortcomings of existing techniques. A description of chosen ice rubble field load events and signal post processing is added. A linear Mohr-Coulomb material model was used for the bulk element. For the cohesive element formulation, a material model was chosen which is based on three irreversible mixed-mode interaction with an arbitrary normalized traction-separation law governed by a load curve. The elastic modulus and fracture toughness for the ice rubble field were scaled using the ice rubble field porosity. A parametric study was conducted, and effects were documented. The numerical model predicted similar values for maximum total force, but average and standard deviation values of total force were higher than measured. The observed load drops in measured force time histories were reproduced with reasonable accuracy in simulated force time histories.


Introduction
The term ice rubble field can be used as generalised term for ice ridge fields where vast area of sea covered with ice blocks packed together. The sea ice ridges are common ice feature in the Arctic and Subarctic sea which are the major load contributor on offshore and marine structures, if icebergs are not present in that area. The ice rubble is an integral part of ice ridges. The ice rubble field is made up of piles of ice blocks and can be fully or partially consolidated. The ice rubble field can be found in between colliding ice floes and can cover large areas. The in-situ testing gives opportunities to understanding the properties and failure of ice rubble in actual environmental conditions. The ice rubble has always been considered as a granular material and can be described by the Mohr-Coulomb criterion, see ISO19906 (2010). Therefore, efforts have been devoted in the past to find reasonable values of the cohesion and the angle of internal friction, see Cornett and Timco (1996), Løset and Sayed (1993), Timco et al. (1992), Sayed et al. (1992), Urroz and Robert (1987), Gale et al. (1987), Heinonen (2004) and Liferov and Bonnemaire (2005). Literature reviews of test methods and analysis of results can be found in Ettema and Urroz (1989), Timco and Cornett (1999), Kulyakhtin and Høyland (2014) and Boroojerdi et al. (2020). But numerical simulation methods have an advantage over laboratory and in-situ testing of ice rubble as they offer a greater number of (design) parameters to be tested under controlled conditions. Therefore, interest in developing numerical simulation methods is increasing. Additionally, simulation methods help to increase the accuracy of the analytical methods used to estimate the design loads.
The numerical simulations enable us to separate and identify the contributions from various parameters controlling the ice rubble field structure interaction process. The ice rubble field is anisotropic material. The ice rubble field presents challenges in terms of constitutive modelling due to, for example, the simultaneous presence of freeze bonds and discrete nature of the ice blocks. The failure behaviour of ice rubble, often called the failure mode, is important information for the estimating of ice loads. The ice rubble failure is dominated by mixed mode failure and the rubble deformation process is controlled by several factors including compressive strength, cohesive strength, fracture toughness and the ice to ice friction. Thus, a robust and reliable numerical method is required to predict realistic load levels created by unconsolidated or partially consolidated ice rubble. The local pressure data can be linked to global ice load, see Bjerkås (2004). However, Fransson and Lundqvist (2006) have proposed a statistical approach to determine the correlations of segmental loads. The failure process is extremely sensitive to the initial conditions. For example, two different failure modes can be observed with and without the ice rubble present at the interaction zone, see Määttänen and Hoikkanen (1990). The failure process of ice rubble is also affected by the parent ice thickness, rubble pile height, ice blocks sizes, the temperature, salinity and the strain rate. Interacting structural properties such as its inclination, width and stiffness, also affects the failure process.
The numerical simulations can be used to study the ice rubble field loads on structures. The challenges involved in the application of numerical methods in ice mechanics are well documented in Bergan et al. (2010). In the simulations of ice rubble field structure interactions, it is necessary to consider a complicated process of fragmentation of ice, formation and motion of ice blocks, and interactions between the blocks as well as between the blocks and the structure. The Discrete Element Method (DEM) has been used to simulation of "fracture prone" and "discrete" material such as ice rubble, see Evgin et al. (1992a), Evgin et al. (1992b), Hopkins (1997) and Paavilainen et al. (2011). In DEM, each ice block is modelled as separate non-continuum element, usually spherical. The forces acting on each element are then computed from the initial properties and the relevant physical laws and contact models. The combined FEM-DEM approach was proposed to eliminate the rigidity issue of bulk element in the simulation of ice rubble structure interactions, see Polojarvi and Tuhkuri (2009), Paavilainen et al. (2011) and Polojarvi and Tuhkuri (2013). In their approach, DEM was used to model the contact interactions and finite elements were used as constitutive relation which dictates the behaviour and the ice fracture. Therefore, it is clear that without solving mesh entanglement problem, conventional FEM cannot be used to simulate ice rubble structure interaction. In contrast, Wong and Brown (2018) have proposed a model of the interaction between ice and conical bridge pier based on FEM where ice sheet is modelled as linear elastic model and considers loads from ice rubble build-up. The recent development in mesh-free formulation techniques of SPH gives an accurate solution for large displacements that remain in continuum domain of Lagrangian framework. For application of SPH formulation to simulate a behaviour of ice rubble in the punch through test, see Patil et al. (2015). The same technique has been used to simulate behaviour of brash ice in an in-situ test by Patil et al. (2020).
A numerical technique called the cohesive zone model (CZM) has been used for analysis of fracture in brittle and ductile materials, see Barenblatt (1962) and Dugdale (1960). The implementation of the CZM into numerical analysis has been termed the Cohesive Element Method (CEM). The CEM requires the insertion of inter-elements between bulk elements in the conventional finite element mesh. Upon reaching the threshold limit stress, the inter-elements (cohesive elements) fail and are removed from the calculation thus enabling the simulation of a "crack" in ice. These newly formed finite element surfaces are then free to interact with each other. The CEM is based on the robust mathematical framework of conventional FEM and is capable of explicitly simulating the fracture process zone (FPZ) which is confined by finite element boundaries. The cohesive element does not have any significant physical mass, thus removing them does not violate any mass conservation law. The cohesive elements describe the cohesive forces in the material when a fracture occurs. The cohesive element method is not new when it comes to simulating interaction between level ice and structure, see Gürtner (2009), Konuk and Yu (2010a), Konuk and Yu (2010b), Daiyan and Sand (2011), Kuutti et al. (2013), Lu et al. (2014), Wang et al. (2018), Wang et al. (2019). The limitations and difficulties associated with the application of Cohesive element method are highlighted by Pang et al. (2015).
In this paper, a numerical model is used to simulate the ice rubble field structure interaction process, using Cohesive Element Method (CEM). The details of constitutive material models used for bulk elements and cohesive elements are given. Both material models are well described and implemented in the LS-DYNA finite element code. A special purpose MATLAB (2019) script was written to post process the load event acquired/measured data, as well to create numerical model and post-process the simulation results. This numerical model is used in two separate parametric studies to investigate the effect of influential parameters such as cohesive strength, angle of internal friction, fracture toughness and exponential decay coefficient. The numerical results are compared with the load event time series and the conclusions are drawn based on merit and drawbacks of the material model and numerical method.

Ice rubble field load event data
Full-scale measurements of the ice loads at Norströmsgrund lighthouse were obtained in two measurement campaigns, "Validation on Low Level Ice Forces on Coastal Structures" (LOLEIF) and "Measurements of Structure in Ice" (STRICE) from 1999 to 2003. The ice forces acting on panels mounted on the lighthouse at water level and environmental variables such as ice thickness, wind speeds, air temperatures, water stage, etc were recorded and documented. Accelerometers, tiltmeters and inclinometers were also used to determine the structural response under ice loading. Additionally, continuous video recording from four different cameras, was done for later interpretation of the data and analysis. A detailed logging of ice condition observations and failure modes was conducted during the measurement periods. The detailed description of the experimental setup at the Norströmsgrund lighthouse is given elsewhere, see Jochmann and Schwarz (2000), Haas and Jochmann (2003) and Bjerkås (2006). Therefore, only a brief introduction will be presented herein. In this paper, two load events were selected from STRICE 2002 measurement campaign data for further study. The selected load events were further post processed to extract and visualise the valuable information.
The Norströmsgrund lighthouse was instrumented with 9 (nine) panels to measure the ice contact forces at the waterline. The panels covered a span of 162 • ; capturing ice loads from the North to the Southeast. The panel no. 9 was oriented to the East and divided into 8 subpanels, which permitted the detection of ice forces at different depths. Subpanels 9-1 and 9-2 make up the top row, followed by subpanels 9-3 and 9-4 in the second row, 9-5 and 9-6 in the third row and 9-7 and 9-8 in the bottom row. In the actual measurement setup, the original panel numbers were not in sequence. To maintain the continuity and to avoid any confusion, in later force analysis new panel numbers (shown in squares) were assigned, see Fig. 1 (a). The ice conditions were monitored as well. The waterline diameter of the lighthouse pier was 7.5 m and the height of base on which the lighthouse is 7.5 m, see Fig. 1 (b). The measured force-time data is used for further analysis and visualization of load events as follows.

Signal processing
Since the load panels were measuring the normal forces only, the resultant force acting on the lighthouse was resolved from the individual total forces measured by the panels. These forces were split into two components referring to a fixed N-E coordinate system of the force measurement system. One of the purposes of this study was to calibrate numerical and material model for these type of load events. Thus, from the simulation point of view, it was better to have a coordinate system which corresponded to the ice drift direction. So, further data processing was done to generate two force components in a X-Y coordinate system, where the Y axis was parallel to the ice drift direction. In this way, ice to structure friction contribution could be added later into force analysis, see Sand and Fransson (2017). This was achieved as follows. The panel forces were split into components referring to a X-Y coordinate system based on ice drift direction angle θ. In this X-Y coordinate system, x-direction was normal to ice drift direction and y-direction was parallel to ice drift direction. The following equations were used to calculate the total resultant force. The vector sum of the load panels normal to the drift direction P x and parallel to the drift direction P y , were defined as: where, n p = 9is the total number of active panels. The contribution from each load panel to the vector sum P x and P y is: where p i (i = 1, 2...n p ) are the measured panel forces acting normal to the face of the panel and μ is the ice to structure friction. Angles α i are measured from the drift direction to centre of each panel as defined in Fig. 1 α i = ϕ + α 0 |n − i|, i = 1, 2..., n p wheren p = 5is the panel number normal to the drift direction, α 0 = 18 • is panel coverage angle and φ is the angle measured from centre of panel n to the drift direction: Angles β i are introduced to keep track on the sign of the contributing forces to the vector sum acting normal to the drift direction: where, n p = 9is the total is number of active panel and n = 5is the panel number normal to drift direction. The total global ice load,P t is calculated from the forces P x and P y acting in x-direction and y-direction.

Sub-load events
Two load events were identified, from the STRICE 2002 database, for further analysis. The load event 1 (identified as 0603_0600 in the measurement database) is from March 6, 2002 database, which starts at 07:24:00 and ends at 07:58:54. The load event 2 (identified as 2103_0900 in the measurement database) is from March 21, 2002, which starts at 20:29:00 and ends at 21:03:00. The vector sum of the load panels in X-Y coordinate system is shown in. Fig. 2 and Fig. 4 for load event 1 (0603_0600) and 2 (2103_0900), respectively. For the simulation purposes, a span of 120 s from force time history was selected for comparison. The sub-load event 1 (0603_0600) is a span of 0-120 s and the sub-load event 2 is a span of 60-180 s. Figs. 3 and 4 gives the summary of ice force distribution on each individual panel for sub-load event 1 and sub-load event 2, respectively. In case of load event 1, despite the ice drift was from 45 • at 0.25 m/s, maximum force was recorded at panel 1 (North). This might be due to pre-existing ice accumulation around the structure. The water stage/level recorded was 0.3 m which means local forces from thinner ice might be missing. This confirms by looking at segmented load panel, see Fig. 3b. In case of load event 2, maximum force was recorded at panel 5 (see Fig. 5a). The load event has avg. ice thickness equals to 0.5 m. The major environmental data recorded for these sub-load events are shown in Table 1. The summary of ice pressure acting on each individual panel for sub-load event 1 (0603_0600) and sub load event 2 (2103_0900) is given in Table 2 and Table 3, respectively.

Constitutive relationships for the ice rubble field
At freezing temperature, varying degree of porosity and freeze bond strength or cohesive strength can exist across the depth of ice rubble field. Often, ice rubble, as a material, has been characterized using properties from soil mechanics such as shear strength, friction angle and cohesion. The shear strength of the ice rubble has four components namely interlocking, frictional contacts, freeze-bonding and failure of the ice blocks, Hopkins (1991). The interlocking and frictional resistance are macro-level phenomena and affected by the shape and the size of the ice blocks. Whereas, the breaking of the freeze bonds and failure of the ice blocks are considered as micro-level phenomena, controlled by tensile, compressive and shear strength. Therefore, both phenomena should be addressed adequately in the material model. A study of Repetto-Llamazares et al. (2011) suggests that a Mohr-Coulomb like failure model can be appropriate for the representation the freeze-bond shear strength as a function of the normal confinement. The Mohr-Coulomb criterion can be used to estimate the strength of ice rubble, see ISO19906 (2010). In the Mohr-Coulomb yield criterion, the ice rubble behaves as an elastic, ideal plastic material. For representative volume for bulk elements, an isotropic elastic-plastic material model was proposed by Hilding et al. (2012) and Hilding et al. (2011). Wang et al. (2018) has reported a decrease in the mean loads when using a linear softening plastic model compared to the perfect plastic model and concluded that contact forces decrease with cohesive softening of bulk elements. However, Kulyakhtin and Høyland (2015) have shown that the Mohr-Coulomb criterion cannot be used to define shear strength of ice rubble where high values of angle of internal friction involved. Despite these limitations, a Mohr-Coulomb yield criterion, linear elastic and ideal plastic constitutive model, can be utilized as bulk elements to  In this study, rubble field geometry was discretized with bulk elements connected to each other by cohesive elements. A material model based on Mohr-Coulomb criteria was used for bulk elements, whereas mixed mode cohesive element formulation was used for cohesive ele-ments, see LS-DYNA (2017). The Mohr-Coulomb criteria was used for solid elements only and intended to represent the ice blocks. This criterion describes the dependence of shear stress τ on the normal stress σ n and is given as:  Where, c and φ are the cohesion and the angle of internal friction, respectively, which controls the fracture stress in mode I and mode II respectively. The compressive C x , tensile T x and shear strength τ are given as follows In the parametric study, the values of cohesion c and angle of internal friction φ were altered to determine their effect on load levels and failure pattern. The Mohr-Coulomb criterion considers effective shear strain as a deformation measurement, which is also in accordance with the assumption of volume preservation.

Cohesive element formulation and constitutive relationships
Granular material such as ice rubble can be effectively analysed for its cracking behaviour by using the cohesive zone method. It involves a gradual reduction of stress or strain in tension or shear loading during crack propagation. In this method parameters describing the crack tip are used to the study of instability phenomena of cracked bodies. These parameters include the strain energy release rate, the fracture toughness, and the crack tip opening displacement (i.e. separation). In the present research, the tension softening is used to analyse the behaviour of ice rubble. For the cohesive element, a material model (*MAT_186) based on three general irreversible mixed-mode interaction cohesive formulations with arbitrary normalized traction-separation law given by a load curve (TSLC), was used to connect bulk elements, see LS-DYNA (2017). This material model was chosen due to the flexibility it provides for the traction separation law. The tabulated traction separation can be defined directly for both fracture modes i.e. I and II, see Fig. 6 (a). The interactions between fracture modes (I and II) are considered and irreversible conditions are enforced through a damage formulation (unloading/reloading path pointing to/from the origin). The Table 1 Sub-load event data from STRICE 2002 database.   traction-separation curve governs the deformation of cohesive elements subjected to tensile or shear stresses. When the force or deformation reaches a critical value, the cohesive element is deactivated or deleted, resulting in crack formation and dissipation or release of the fracture energy G F . A crack can grow via deformation and failure of neighbouring cohesive elements. The traction separation behaviour of this model is mainly controlled by the energy release rate and peak traction. The energy release rateG c I and peak traction T in normal direction defines mode I, while, the energy release rate G c II and peak traction τ in tangential direction defines mode II. The maximum separation δ F I and δ F II in mode I and mode II are given respectively as: whereA TSLC is the total area of the normalized traction-separation curve, as shown in Fig. 6 (b). In the current research, for mixed mode formulations of cohesive element, the effective separation parameter TES and exponent of the mixed mode criteria XMU were set to 0 and 1, respectively. The ultimate mixed-mode displacement δ F (i.e. total failure) for this formulation of cohesive element is given as: where β = δ II /δ I is the ratio of mode mixity. In this model, damage of the interface was considered, i.e. irreversible conditions are enforced with loading or unloading paths coming from or going to the origin. The total mixed-mode relative displacement δ m is defined as where δ I = δ 3 is the separation in normal direction (i.e. mode I) and is the separation in tangential direction (i.e. mode II).
To determine the shape of the traction separation curve (TSLC), a three-element setup was used. The detailed description of calibration procedure was given in Sand and Fransson (2017). In this setup, two bulk elements are joined together by a cohesive element. Then, the bulk elements are subjected to uniaxial tensile loading (i.e. pulling away from each other). This causes the tensile softening of cohesive element as crack is assumed to occur between two elements. The traction-separation curve as shown in Fig. 6 (b) is based on damage formulation that accounts for strain softening. The damaged traction t d i in mode I vary with the nominal tractiont i in i-direction(i = x, y, z) as follows: The damage parameter d i increases from an initial value of zero to a maximum value of d max = 0.99. Moreover, damage in the cohesive elements accumulates with traction t i in mode I, in accordance with the following function: where T i is the threshold value for damage in mode I in i thdirection(i = x, y, z) and D is a softening parameter set to 1.0. The mesh size sensitivity of these elements was controlled by maintaining constant fracture energy, regardless of the element size. This was achieved by including the equivalent element length 'L eqv ' (cube root of the hexahedral element volume). Similarly, the shape of the traction separation curve in mode II can be obtained by assuming the same fracture toughness as in mode I (i.e.K C I = K C II ) and threshold fracture stress τ. The effective elastic modulus E eff for cohesive elements and is set equal to5/6E. In this way, a cohesive element stiffness will always be lower than for a bulk element. The cohesive elements are act as a potential crack plane and deformation occurs according to a given tractionseparation curve, when cohesive elements are subjected to tensile or shear stresses. Higher cohesive element stiffness can lead to numerical instability in explicit time stepping schemes. In the mixed mode fracture where loading is both in normal and tangential directions, a failure envelope type criterion that considers all stress and damage components can be utilized for more accurate fracture predictions.

Estimation and scaling of material properties
In the rubble field structure interaction, the value rubble porosity affects the magnitude of rubble load. The rubble porosity η r is not measured, and value varies between 0.2 and 0.5 in various literature including Høyland (2000). Assuming value of porosity close to 0.5 will underestimate the load as much of energy is absorbed by rubble itself. Whereas, value close to 0.2 will results in higher load due to local crushing and shearing at the interaction zone. For this analyses purpose ice rubble density ρ r is assumed to be 700 kg/m 3 which gives porosity value of 0.24. The mechanical properties of ice rubble field are strongly dependent on the properties of the parent level ice sheet. Therefore, it is likely to develop a relationship between properties of parent level ice sheet and ice rubble field. Fransson and Stehn (1993), has shown the effect of porosity η r on the elastic modulus and the fracture toughness of the warm ice. Using same the analogy, as a preliminary approach, properties of parent ice sheet were scaled by factor of (1 − ̅̅̅̅ η r √ ) to obtain properties of the ice rubble field. The stiffness depends on the pressure rate in granular materials with high porosity such as the ice rubble field. In this study, following scaling was used to estimate the elastic modulus of the ice rubble field E r based on elastic modulus of ice E ice .
Mulmule and Dempsey (1997) have studied mode I fracture of sea ice by using a fictitious crack model. The behaviour of TSLC is mainly controlled by fracture energy (G c I for mode I failure and G c II for mode II failure), fracture stress (T for mode I failure and τ for mode II failure) and the shape of TSLC. The fracture energies (G c I &G c II ) are defined as the energy release rate to create a unit crack and can be obtained from fracture mechanics tests for level ice. Based on in-situ tests, Dempsey et al. (1999) estimated fracture toughness for ice as between 0.170 MPa ̅̅̅̅ m √ for small ice specimen and 0.25MPa ̅̅̅̅ m √ for large ice floes. Furthermore, according to Dempsey et al. (2004) and Dempsey et al. (2012), the size-independent fracture energy of first year ice is 15 J/m2, and for multi-year (MY) ice the value is in the range of 23-47 J/m2. However, the uncertainty related to scale effect of test specimen is reported by Timco and Weeks (2010). In the simulations, the fracture toughness varied between 0.100 and 0.250 MPa ̅̅̅̅ m √ . It should also be noted that by using a relatively low value of fracture toughness the numerical model will become too brittle and too weak.
The critical fracture energy of the ice rubble field, G F r , can be expressed in terms of fracture toughness, K c Ir , poison's ratio, ν and elastic modulus E r as follows Similarly, scaling of the fracture toughness for the ice rubble field based on its porosity can be done as follows.
So, by using Eqs. (14) and (16) the fracture energy in Eq. (17) can be rewritten as follows As the fragmentation of rubble field depends on complete dissipation of the fracture energy, the fracture energy will be one important parameter that needs to be estimated as accurately as possible. In the case of the rubble field, the yield strength is very low compared to the level ice, which means the failure of the ice rubble field mainly depends on the cohesive element properties. A parametric study on the cohesive strength and fracture energy was performed to investigate the effects on load level, failure mode and behaviour of ice rubble field. The proposed scaling for properties of ice rubble field can be basis for further investigations.
In CEM, the TSLC is used to define the relationship between traction (T for mode I and τ for mode II) and separation δ, as given in Eq. (8). The failure type (ductile or brittle) is mainly dependent on the shape of TSLC. Two factors affect the shape of TSL, namely, the fracture stress and fracture energy. The fracture stress is defined as the threshold stress to initiate a crack. The three most common types of cohesive laws, linear, exponential and trapezoidal softening, are commonly used in engineering practices. To model the ductile type fracture, a trapezoidal softening curve is used. This type of curve has a softening behaviour with a significant plateau region after the peak. To model brittle crack growth, linear and exponential models are often used. They give less fracture energy than the trapezoidal softening law, as the area under the TSLC represents the fracture energy per unit area, i.e. the amount of energy dissipated when the crack was fully developed. However, based on experimental investigation of elastic-plastic solids such as aluminium, Tvergaard and Hutchinson (1992) and Cornec et al. (2003) have concluded that the effect of the TSLC shape is insensitive to the resulting fracture behaviour as long as the fracture energy of the model is correct. There are not any fixed and definitive criteria to select the shape of TSLC as the CEM is only an approximate mathematical representation of the physical fracture process. Modelling of the fracture process of ice rubble field is complicated due to high porosity and presence of cohesive bonds. The TSLC shape is affected by the change in fracture toughness of level iceK c I , shear strength τ and angle of internal friction φ. The TSLC used in the parametric study associated with sub-load event 2 are given in Fig. 7 to Fig. 9. The TSLC with three different values of level ice fracture toughness used in scaling for mode I is given in Fig. 7. The TSLC with two different values of fracture stress (i. e. T i ) for mode I is given in Fig. 8, whereas, the TSLC for mode II with three different values of fracture stress (i.e. τ), are given in Fig. 9.

Numerical modelling of ice rubble field interacting with the Norströmsgrund lighthouse
The finite element method with Lagrangian element formulation was chosen to simulate the lighthouse ice rubble field interaction load event. The commercially available explicit finite element software LS-DYNA, version R10, was used for modelling of ice rubble field and lighthouse. The finite element model used to discretise the rubble field geometry had two different parts; cohesive and non-cohesive, as shown in Fig. 10 (b). The cohesive elements were only used to discretise the middle part of ice rubble field, termed as the cohesive part. The cohesive part was supported by bulk elements with coarser mesh densities, termed as noncohesive part. This arrangement is useful in saving computational time as the number of nodes in the non-action area are reduced. Also, the purpose of the non-cohesive part was to provide a confinement for the cohesive part. The cohesive part was discretized with standard solid finite element e.g. hexahedral elements. Each individual solid element was connected to its neighbouring solid element using cohesive elements. As finite thickness cohesive elements were found to perform better numerically than zero-thickness elements, see Kuutti et al. (2013), a very small thickness (=10 − 4 mm) is given to cohesive elements.
The shear strength of ice rubble is highly dependent on the boundary conditions, see Fransson and Sandkvist (1985), Ettema and Urroz (1989), Hopkins (1991) and Timco and Cornett (1999). Therefore, realistic boundary conditions are necessary to simulate load levels accurately. The dimensions of ice rubble field used in simulations are shown in Fig. 10 (b) and are given in multiples of the diameter of structure (D = 7.5 m) at water line. The dimensions in X and Y directions were kept five and six times the diameter of structure, respectively. The global axis system is shown in Fig. 10. The dimension in Z direction was equal to the average thickness of the load event. The outer nodes of ice rubble field were fixed in Y and Z direction. The nodes at the edge of the opposite side of structure and side edges of ice rubble field geometry were fixed in X and Y directions, respectively. The preliminary studies confirmed that the chosen dimensions of ice rubble field were not affecting the interaction process and load levels. The lighthouse structure was modelled with 1:1 ratio using shell elements as per dimensions from STRICE 2002 database. The bottom nodes at foundation of the lighthouse structures were fixed in all directions. The general geometry and simplified finite element model used in the parametric analysis are shown in Fig. 10. The basic properties associated with ice rubble field are given in Table 4. The elastic material model is used for lighthouse parts and the material properties used are given in Table 5. Furthermore, ice rubble field was moved toward the structure with constant drift speed of the load event. The model generation process was automated by using a special purpose MATLAB script. The following assumptions were made in the finite element model.
1. Temperature variation in ice rubble is neglected i.e. ice rubble is assumed to have constant temperature. 2. No spatial variation is considered in ice rubble i.e. ice rubble is assumed to have same material properties throughout.
3. The volume of the ice element is conservative during the deformation. 4. The strength of the rubble is uniform which means there is no spatial variation in the strength of the rubble field. 5. Thickness of the ice rubble field is equal to the average ice thickness of the load event.
The fully integrated solid hexahedral element formulation was used for a bulk element, while, 8-node cohesive element formulation (ELFORM = 19) was chosen for a cohesive element. The main drawback of hexahedral element type of mesh is "zig-zag" crack pattern, see  Konuk and Yu (2010b). In this type of mesh, regardless of any mesh refinement, a crack travel path is about √2 times longer than the intended crack path, leading to an extra energy consumption. The total error should be close to 1 − ̅̅̅ 2 √ . The error can be reduced by using triangular-shaped or tetrahedron mesh but with an increased computational time. Also, it was reported by Wang et al. (2018) that the mean load increased with mesh refinement and convergence was not achieved. In order to obtain realistic global horizontal loads, Lu et al. (2014) suggested that, mesh size preferably be twice the ice sheet thickness and has also shown that the fragmentation size is largely dependent on the element size. A sufficiently small mesh size is required to simulate microcracks. However, computational time increases rapidly with smaller mesh size. So, it is    A. Patil et al. impractical to use very small mesh size with cohesive elements to include the microcrack effect. To address this issue some other techniques must be implemented for example an adaptive mesh refinement. The crack pattern in ice rubble field and structure interaction is somewhat unknown and to keep the computational cost down, hexahedral elements were used. The mesh size was kept constant i.e. mesh sensitivity was not part of this study. In this paper, bulk element size for load event 1 (0603_0600) was 200, 200 and 190 in X, Y and Z direction respectively and for load event 2 (2103_0900), bulk element size was 200, 200 and 125 in X, Y and Z direction respectively. In the numerical analysis, the contact forces acting on the load panels were extracted by using force transducers. Note that the contact forces on the load panels were written in the global coordinate system and were transformed in direction normal and tangential to the load panels.

Hydrodynamic forces
The buoyancy and drag forces of water are very important components in ice structure interaction simulation. There are several methods available to include the buoyancy and drag into the numerical model, for example the arbitrary Lagrangian-Eulerian (ALE) and CFD, in which the water is represented as an explicit body in simulations. But this explicit representation is very time consuming hence costly. In the present study, a cost effective way to include the buoyancy and the drag forces of water, is presented by using a discrete mass-spring-dashpot system, shown Fig. 11. The same method is used by Sand and Fransson (2017) to simulate level ice structure interaction with a lighthouse using cohesive element method. A spring-dashpot is connected to each node of the eight nodes of hexahedral element. The buoyancy force F b,i is a function of the displacement Z t relative to waterline as shown in Fig. 11(b) and (c). The buoyancy force on each fully submerged node of hexahedral bulk element of length L ex , L ey and L ez in X, Y and Z direction respectively, was calculated as follows: where, F b,i , ρ W and gare buoyancy force, density of the water and gravitational acceleration, L x , L y and L z is the length of the hexahedral element in X, Y and Z-, direction, respectively. The drag force F d,i is defined by a basic viscous damping equation for an object moving with a vertical velocity V through a liquid: where C D is the drag coefficient and in all simulations the value of 1.05   was used which is equal to a drag coefficient of cube moving through a fluid. Then, the total forceF T,i on each node for the discrete mass-springdashpot system in global Z-direction is given as where, ρ r is density of ice rubble.

Frictional forces
When the force or deformation reaches a limit value the cohesive element is deactivated (deleted) from calculation. Once they are deleted, free solid elements may interact with each other with pre-defined contact condition. In LS-DYNA, the frictional coefficient μ is assumed to be dependent on the relative velocity V rel of the surfaces in contact and calculated as follows where, μ S , μ D and Dc are the static, dynamic and exponential decay coefficient of friction, respectively. The friction of saline ice was measured in a laboratory by Kennedy et al. (2000) and Frederking and Barker (2002) and in the field by Sukhorukov and Løset (2013). The laboratory measurements by Kennedy et al. (2000) were done at low sliding velocities and with both saline and freshwater ice. However, the results of those experiments were for low sliding velocities only and do not indicate significant differences between the friction coefficients of saline and freshwater ice. Fortt and Schulson (2007) suggested that peaked shape of friction curves versus sliding velocity was due to a change in sliding behaviour. Sukhorukov and Løset (2013) studied the friction of sea ice on sea ice performed in the Barents Sea and fjords at Spitsbergen. The presence of sea water in the sliding interface had very little effect on static and kinetic friction coefficients. The measurement conducted in the field by Sukhorukov and Løset (2013) give much higher friction coefficients than those obtained in the laboratory by Kennedy et al. (2000). This might be due to the fact that in the field the environment conditions are not easy to maintain, and there might be additional processes like ploughing involved due to a higher roughness of the ice surface. Fransson et al. (2011) showed that surface roughness, speed and temperature, play an important role in deciding the friction coefficient on ice surfaces. Similarly, Tikanmäki and Sainio (2020) concluded that the surface roughness is the most influential parameter regarding the friction coefficient based on laboratory experiments. The ice temperature was not measured for load events. Therefore, an average of air temperature and water temperature values was used to estimate the ice temperature of ice for the sub-load event. For the sub-load event 1 (0603_0600) and for sub-load event 2 (2103_0900), the ice temperature T r was estimated to − 6 • C and − 0.65 • C, respectively. Based on estimated ice temperature, the decay coefficients were selected in such way that, a fit was obtained to reported values of test data. The static, dynamic and decay friction coefficients used in simulations are given in Table 6. Modelled ice-to-ice friction coefficient vs values reported by others are shown in Fig. 12.
Two types of penalty contact formulations were defined. An eroding single surface contact formulation was used for ice block to ice block contact due to failure of the ice sheet. The ice block-to-ice block friction model mimics sea ice to sea ice friction, as described in Eq. (21). For ice block to lighthouse contact, an eroding surface-to-surface contact was used. For this contact, the static friction μ S , the dynamic μ D , and the decay D c coefficient were set to 0.1, 0.05 and 0.02 respectively.

Damping
When modelling ice structure interaction with CEM the contact stiffness in the contact definition cannot be ignored, see Feng et al. (2016). The kinetic energy of impacting ice sheet is converted into ice fragmentation which lead to bulk elements being expelled at high velocity. This is a non-physical phenomenon and results in predicting higher forces. To obtain a stable model, the viscous damping coefficient VDC was used in contact definition. The VDC damps out unwanted oscillations normal to the contact segment. The VDC should be an integer and expressed in percentage of critical damping. The value of 20% (recommended by LS-DYNA for impact simulation) was applied in all performed simulations.
The lighthouse structure was modelled with the ratio of 1:1 and used material properties for steel, concrete and sand. Consequently, an appropriate damping force is needed to mimic the real-world structural response to external forces such as collision with ice. For the numerical lighthouse structure, the damping frequency range deform method was chosen. The user must input three parameters to define this damping. For all numerical simulations, 2% of the critical damping and the lowest (f low ) and highest (f high )frequencies were set to 2.5 Hz and 200 Hz respectively. Also, the mass-weighted damping was used to damp out the unrealistic motion of fragmented bulk elements. As a preliminary approach, the lowest eigenfrequency of the structure (f low = 2.5Hz) was used in calculation of the damping coefficient as follows The damping ratio ζ was chosen to be 0.75% and 0.5% for sub-load event 1 (0603_0600) and 2 (2103_0900), respectively.

Analysis of numerical simulation results
The series of numerical simulations in the parametric study were conducted to identify and to quantify the influence of the parameters on simulation results. The most important aspects of simulation results were the deformation of the ice rubble field, failure pattern (i.e. crack pattern), the transportation of ice pieces from the crushed zone, ice rubble formation process, the piling of ice in front of the lighthouse, the simulated force level and nature of force. Two parametric studies were performed, one for each load event, showing the effects of cohesion c, angle of internal friction φ, fracture toughness and friction (decay coefficient in particular). The measured force time history of two full scale sub-load events were simulated using the CEM. The simulation time span was 60 s and was sufficient to capture the overall trend. All the computations are done using double precision solver provide by LS-DYNA© R10.0 with 8 cores (CPUs) running in parallel on a server and supervised by a real-time monitor. The explicit time integration scheme was applied. The mass scaling was applied to increase the stable time increment. No significant increase of the mass of the model was observed due to the small thickness of the cohesive elements. Bulk elements experienced plasticity while being compacted in the failure region. For sub-load events 1 and 2 the data were sampled at frequency of 50 Hz and 30 Hz, respectively. The sampling frequency for all simulations was kept at 50 Hz.
At the start of simulation (t = 0) the recorded force was zero, as there was no contact between ice and lighthouse. This differs from reality as a measured load event may begin with extra load on lighthouse due to unknown initial condition. Also, the reported water level for load event 1 (0603_0600) was above the load panels and initially ice was hitting only part of panels, which makes the ice load highly sensitive to initial conditions Consequently, it is very difficult to mimic the initial conditions in simulations. In simulation, the initial condition is always zero loads or no contact, while the measured load events always begin with load that has already built up on structure. Despite this limitation, the simulation 1-1 was able to capture the dynamic nature for forces as well as ice accumulation behaviour. The load event 2 has smaller average ice thickness than load event 1 (0603_0600) which lead to lower panel loads. In simulation 2-1, estimated material properties (i.e. cohesion, angle of internal friction and fracture toughness) has yielded load peaks that match to the measured ones, see Fig. 14.
As the simulation progressed, more and more ice came in contact with lighthouse, which caused an increase in the force, eventually reaching a stable level. The load panels no. 5 and 6 were showing highest force levels as they were aligned to the drift direction, see Fig. 17 and Fig. 18. The frequent force drops can be seen in total panel forcetime figures (Figs. 13, Figs. 14 and 19-25). The hexahedral elements, used for bulk elements, acted as small cantilever beams due to finite length and stiffness. As deformation continued these small cantilever beams bended and buckled due to contact and friction. The compressive crushing process was a cyclic process which causes these load drops. The ice edge at the contact surface was moving either upwards or downwards during the loading and unloading phases. The repeated vertical motion of the ice edge was caused by an imbalance in the counteracting forces that arise from the rubble pile formed above and below the ice sheet. The ice rubble tends to form under the ice sheet, near the crushed zone until it clears out.
Two separate parametric studies were conducted for each sub-load event. In each simulation only one parameter was changed at a time thus enabling identification of the effect of that parameter both qualitatively and quantitatively. Input parameters of parametric studies for simulating the sub-load event 1 (0603_0600) and 2 (2103_0900) are given in Table 7 and Table 8, respectively. The comparison between measured and simulated global forces P t for sub-load event 1 (0603_0600) and 2 (2103_0900) are presented in Figs. 13 and 14, respectively. Partially enlarged screenshots of the ice rubble field deformation in the front of the lighthouse for simulations 1-1 and 2-1 at t = 2, 12, 20, 37, 50 and 60 s are shown in Figs. 15 and 16, respectively. Both simulations 1-1 and 2-1, captured the dynamic nature of the measured forces as well as the ice rubble formation process. The detailed input to simulation 1-1 and 2-1 is given in Table 9. The force drops near the end of simulation 1-1 (see Fig. 13) which could be linked with clearing of rubble in front of lighthouse foundation. In the simulation of sub-load event 1 (0603_0600), the lighthouse penetrated 15 m into the ice sheet at the end of the numerical analysis. In simulation 2-1, the Fig. 12. Friction law implemented in LS-DYNA and fitted to data for sea ice to sea ice friction. Fig. 13. Comparison between the total panel forces of the sub-load event 1 (0603_0600) with the numerical simulation 1-1 (c = 67kPa, φ = 40 • , K c force drops were not as high as simulation 1-1, and the relatively narrow scatter of force was registered. The occurrence of peaks and drops of the force can be associated with compressive failure, contact friction and clearing of rubble in front of load panels. Initially, the edge of the ice rubble field hits the structure, which developed a local stress zone thus causing the ice to fail by crushing failure. As the structure penetrated the ice rubble field, ice load built up on several load panels and loose rubble started to form in front of the structure. Once the cohesive elements are deleted, which further leads to further crushing failure, the process of ice rubble formation starts. As the simulation progressed, ice rubble accumulated above and below water line. The rate of ice rubble accumulation was different in the simulations of the two sub-load events due to difference in ice drift speeds. As a result, at the end of 60 s the simulation 1-1 produces more ice rubble  than simulation 2-1. The formation of crushing failure is a cyclic process and causes a cyclic loading response. The continuous crushing failure of the ice sheet and clearing of ice rubble leads to the drops and rises of forces. The vertical and horizontal motion of ice rubble was influenced by dynamic forces from friction and the ice rubble pile above and below the ice sheet. The measured individual load event panel force data was also compared with numerical simulation 1-1 and 2-1 in Figs. 17 and 18  respectively. At the start of the numerical analysis, panel pressure normal to panels 5 and 6 increased rapidly. This pressure fluctuation corresponds to the repetitive occurrence of the crushing failure. The panel 1 was nearly inactive during the simulated time. The measured panel load distribution was somewhat different from the performed numerical simulations. The maximum panel pressure occurred on the panels facing the drift direction i.e. panel 5 and panel 6, while the panel pressure was considerably lower on the panels at the outer edge. In case of simulation of load event 1 (0603_0600), the actual measured load applied on the panel no. 1 was significantly higher than the simulated. It can be related to the unknown initial condition of the load event.

Effect of change in cohesion and internal angle of friction
In the present material model for cohesive elements, the threshold fracture stress, in mode I and II, are based on relationships given in Eq.
(7) which are functions of cohesion c and angle of internal friction φ. It is assumed that the ice sheet representing the ice rubble field is elastic before reaching the threshold fracture stress. After the first crack, the ice rubble field is more prone to develop cracks in the vicinity of the initial crack. Once the crack is complete i.e. fracture energy is reached, the ice rubble field behaves more like a viscous fluid. When threshold fracture stress (in mode I or mode II) decreases, keeping same fracture energy, the traction separation curve tends to exhibit ductile type fracture. This results in narrow scatter of forces. In the case of the decrease in the fracture stress in mode II (i.e. shear strength), "cracks" go sideways as opposed to the case of decrease in fracture stress in mode I (i.e. tensile strength). The values used for cohesion and angle of internal friction were arbitrary and chosen to fit the peak load levels from load event. A comparison of total panel forces of load event 1 (0603_0600) and  Fig. 19. Comparison of total panel forces of load event 1 (0603_0600) with numerical simulation 1-1 and 1-2.
A. Patil et al. numerical simulation 1-1 and 1-2 is shown in Fig. 19. The decrease in cohesion from 67 kPa to 36 kPa (46% decrease) resulted in a decrease of total panel forces by 27% for sub-load event 1 (0603_0600), Fig. 19. Whereas, the comparison of total panel forces of load event 2 (2103-0900) to numerical simulation 2-1 and 2-2 is shown in Fig. 20. While, increase in cohesion from 36 kPa to 67 kPa (86% increase), mean forces increase by 124% as shown in Fig. 20, for sub-load event 2 (2103_0900). The angle of internal friction φ controls the tensile strength T and thus the fracture stress in mode I. As the angle of internal friction φ is inversely proportional to tensile strength, its decrease results in an increase of the tensile strength. The effect of a change in the angle of internal friction on the total panel force for the sub-load event 1 (0603_0600) and 2 is shown in Figs. 21 and 22. Three input values of 20 • , 40 • and 60 • were used. In simulation of sub-load events 1 and 2, a wider scatter of forces is registered with increase in angle values. In comparison, the simulation 1-4, registered highest peak forces, shown in Fig. 21 and the simulation 2-4 in Fig. 22, registered the highest peak forces. Whereas, the lowest value of angle of internal friction (i.e. 20 • ) registered a narrow scatter of total forces.

Effect of change in fracture toughness
The effects of changes in the fracture toughness on the numerical results compared with the load event were investigated. The scaling of fracture toughness of level ice was done by using scaling formula presented in Eq. (16). The range for fracture toughness of level ice used was from 0.100 to 0.250MPa ̅̅̅̅ m √ . Then the fracture energy was estimated based on Eq. (17). The comparison of different values of fracture toughness for level ice used in scaling in terms of total force time histories is shown in Figs. 23 and 24. The fracture toughness has a significant effect on the TSLC. The increase in fracture toughness, keeping elastic modulus same, the total fracture energy is also increased, see Eq.

Effect of change in decay coefficient
When cohesive elements are deleted, bulk elements become separated from the parent ice sheet. Then the contact condition can be applied to implement the frictional forces. The frictional coefficient was assumed to be dependent on the relative velocity V rel of the surfaces in contact, see section 4.3. The input values used in parametric study for static μ S , dynamic μ D and decay Dc coefficients are shown in Table 6.
Two sets of input values of exponential decay coefficients are used in parametric study number. 2, associated with load event 2 (2103_0900). By changing the decay coefficient D c , the range of dynamic friction coefficient can be changed. The values of decay coefficients used in set 1 and set 2 equals to 6.61 and 2. The effect of change of the decay coefficient was only investigated for sub-load event 2, while for sub-load event 1 (0603_0600), the value was kept constant at 10.4. Fig. 25 shows the effect of two different values of decay coefficient on the modelled total force time history for the sub-load event 2 (2103_0900). In addition to the shortcomings of the hexahedral element mesh, described in section 4, this type of mesh has also a larger contact surface area leading to a higher frictional force component. In comparison with the measured sum of panel forces plotted in Fig. 25, the simulation 2-5 yielded higher dynamic friction coefficient, which resulted in higher frictional forces.

Panel line load
The comparison of the spatial distributions of panel load for the load event 1 (0603_0600) and 2 (2103_0900) with numerical simulations is  A. Patil et al. shown in Figs. 17 and 18, respectively. The peak load, mean load and standard deviation for each simulation in the parametric study is presented in Table 10 and Table 11. The distribution of average, maximum and standard deviation of the panel load values for load event 1     Fig. 26 and Fig. 27. Due to higher sampling rate, despite the fact that the individual registered peak lasted for very short time, it can give much higher load than the measured one.

Summary and conclusions
The cohesive element method (CEM) was used to simulate the interaction between the ice rubble field and Norströmsgrund lighthouse. An explicit scheme was used to solve the finite element model in LS-DYNA. The ice sheet was modelled using hexahedral elements (bulk elements) with the Mohr-Coulomb (linear elastic-ideal plastic) material model and the cohesive elements were inserted in between the bulk element mesh by duplicating nodes along all internal mesh boundaries. The measurement data of selected sub-load events was post processed to enable their comparison with simulation results. A parametric study of simulations was performed and the effects of change in various parameters were documented. The identification of the ice rubble structure load event was crucial. The load events were chosen to be simulated based on logbook description, drift angle and fairly constant ice thickness. The drift angle was 45 • (with respect to North direction) for both load events. The measured average ice thickness and the drift speed were 0.760 m & 0.25 m/s for load event 1 (0603_0600) and 0.5 m & 0.15 m/s for load event 2 (2103_0900).
The objectives of this study were to model the ice rubble field lighthouse interaction with cohesive element method and to study the effects of the various influential parameters on simulation results. The numerical model used in the parametric study consists of hexahedral bulk element mesh connected by the cohesive elements and beam elements connected to each node of ice rubble field for buoyancy and drag. The procedure of calibrating material models for bulk and cohesive elements require extensive sets of experimental data for example compressive strength tests, shear strength tests and fracture toughness tests. The absence of such experimental data requires the reliance on assumptions. Thus, scaling of level ice material properties such as elastic modulus and fracture toughness was done to find suitable material properties for the ice rubble field. The chosen scaling yielded reasonable material properties for the ice rubble field. The average ice sheet thickness was estimated from laser and EM measurements. Therefore, the chosen ice thickness did not represent the natural variations in the ice sheet thickness correctly which may affect the accuracy of the model. The assumption of homogeneity of material properties of ice rubble field may not be valid as ice sheet properties vary from within the ice rubble field.
Based on the comparison of the time series for the global loads the influence of simulation parameters was shown. However, the presented numerical model has certain shortcomings regarding simulation of the nature of forces, panel load levels and the fracture behaviour of the ice rubble field. The general conclusions drawn here are based only on the cases studied in this paper and given below.
• The discrete mass-spring-dashpot model is proposed to simulate the buoyancy and drag of ice blocks in water. However, accuracy of the discrete mass-spring-dashpot model is the subject of discussion. The use of cohesive elements only for cohesive zone part of numerical  model can reduce computational time without losing simulation accuracy. • The damping of the lighthouse was necessary to give realistic structural response under ice loading. In the present study, the damping method was applied with an ad hoc approach. Therefore, the influence of the damping method on the failure mode of ice rubble field is not investigated here. The damping effect of water on lighthouse is not modelled. Thus, its contribution towards failure modes and load levels of ice rubble field is unknown. • The mixed mode formulations for cohesive element is successfully implement in this model to simulate ice rubble field structure  interaction. Also, the two-element setup for calibration of cohesive elements is proposed. • The scaling used to estimate the ice rubble field properties, was based on linear scaling factor of (1 − ̅̅̅̅ η r √ ). The ice rubble field can have different porosity profile along the depth. Therefore, the depth dependent material properties cannot be scaled in this way, especially that the dependency between the fracture toughness and porosity is unknown. Consequently, more research is needed for finding an appropriate scaling factor. • The used bulk element material model (i.e. Mohr-Coulomb criteria) did not have any softening feature implemented, which may lead to an overestimation of the predicted forces. • The presented numerical model was able to capture the features of force time history such as load drops and peaks. The rubble accumulation around the structure was also reasonable. However, measured data was not available to compare the rubble pile-up dimensions. • The finite element size used in simulations is small enough to capture the macro-cracking. However, to capture the micro-cracking phenomena, finer mesh size would be required. Furthermore, hexahedral mesh tends to predict a higher load indicating the need of mesh shape and size reconsideration. Thus, a tetrahedral element mesh can be adopted for future research to simulate realistic crack pattern. • The cohesion (c) value has direct effect on load levels. Higher cohesion values tend to give higher load levels. Whereas, change in angle of internal friction has opposite trend. The lower values of angle of internal friction (ϕ) produced higher load levels. The nature force scatter is also changing with higher values of cohesion (c) giving wide scatter. In contrast, narrow force scatter was observed with higher values of angle of internal friction (ϕ). • No general trend was observed with change of fracture toughness. In case of load event 1, higher values of fracture toughness results in narrow scatter of force. Whereas, in case of load event 2, higher values of fracture toughness registered highest forces with wider scatter of forces. • The decay coefficient used in set 1 (D c = 6.61) yielded lower dynamic friction coefficient which registered lower forces compared to set 2 (D c = 2).
However, it should be mentioned that the presented model is at a preliminary stage. Further developments, as already suggested above, are needed. The simulation results showed potential to use this method to simulate ice interaction with other types of structures such as ships or marine vessels. Also, this is a promising first step in showing the potential to simulate other types ice features for example ice ridgestructure interaction. The efforts to simulate various other ice load events will be furthermore be pursued.