Armelle Decaulne

Najat Bhiry

Janie Faucher-Roy

Clara Pelletier Boily

The development of Kangiqsualujjuaq and the threat of snow avalanches in a permafrost degradation context, Nunavik, Canada

Keywords: Nordic village, slopes, active layer, demographic growth, snow avalanches, Canada village nordique, versants, couche active, croissance démographique, avalanches, Canada

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

Introduction 1

In Nunavik, few studies have focused on slope processes (snow avalanches, landslides, rockslides, etc.). Indeed, according to Germain [2016], Nunavik has favourable environmental conditions for gravity movements : steep, bare slopes heavily weathered by strong periglacial conditions, snowfall and strong winds. Recent studies have shown that gravitational processes began in Nunavik about 4600 years ago due to the release of pressure accumulated in rock material during the glaciation period [START_REF] Decaulne | Geomorphic evidence of slope dynamics on the Canadian Shield throughout the Holocene -a study from Lac à l'Eau-Claire, western Nunavik[END_REF], Bhiry et al., 2019]. In the Umiujaq region, southwest of Nunavik, gravitational movements pose a significant risk to the community on the road to Lake Tasiujaq, which passes below the steep high slopes of the Hudsonian cuestas [Veilleux et al., 2019]. Several snow-avalanche events occurred in the past, in Inukjuak, Ivujivik, Kangirsuk, Kangiqsujuaq, Quaqtaq and Salluit [START_REF] Lied | Avalanche hazard assessment in Nunavik and on Côte-Nord, Québec, Canada[END_REF]. However, the deadly snow avalanche that occurred on the night of December 31 st , 1998-January 1 st , 1999 in Kangiqsualujjuaq, located in north-eastern Nunavik (Figure 1) [Bérubé, 2000], highlighted the level of risk. At this time, the inhabitants were gathered to celebrate New Year's Eve in the school gymnasium. The gymnasium was located within the deposit zone of a short snow-avalanche path. The memory of this event is locally longlasting. However the accurate perception of hazard is impeded by the lack of systematic data collection regarding slope activity in locations where hazard could easily shift to risk due to the vulnerability of settlements or short transportation corridors around settlements. In the 1980s and 1990s, at least four snow avalanches occurred in Kangiqsualujjuaq before the 1999 dreadful snow avalanche. Following the 1999 catastrophic snow avalanche, the spatial organization of the village changed significantly, starting with the relocation of several buildings (school, church, and a Co-op store) and houses. The erection of a 150-m long and 5-m high protection wall at the base of the slope was carried out on decision of the Ministry of Public Security, and a 100-m exclusion zone was established between the foot of the slope and the village [START_REF] Schaerer | Avalanche hazard mitigation for Kangiqsualujjuaq, Québec[END_REF]. The subsequent development of Kangiqsualujjuaq in response to its population growth must take into account the permafrost degradation on the land on which it is built. The main objective of this paper is to examine how a single snow-avalanche event caused the breakup of the original structure of Kangiqsualujjuaq village, far from its historical centre, in a context of demographic growth. To achieve this goal we conducted a literature review, a diachronic analysis of aerial photographs and our own geomorphologic observations in the field.

Prone snow-avalanche weather conditions in Kangiqsualujjuaq valley 3 Snow avalanches are a natural phenomenon resulting from the combination of slopes and snow, the latter being destabilized and moving downslope under different conditions. The hazard is well known in Canada, where several recreation, transportation, property and resource industries are at threat [Stethem et al., 2003]. In Nunavik, several villages and village surroundings are at risk, although the issue lacks in-depth studies, as several of the communities gather favourable terrains, i.e., slopes with inclination over 10-15° [Lied and Doomas, 2000], coupled with favourable weather conditions (temperature, precipitation and wind). The issue is exacerbated by weather conditions conducive to the occurrence of avalanches such as strong winds, snowfall and bare slopes heavily weathered by heavy periglacial activity ; however scarce and only recent weather data are available due to several newly installed and partially operating weather stations. The closest precipitation and air temperature measurements from Kuujjuaq indicate the clear polar climate (Figure 2A), with seven months recording a negative mean temperature (period 2006-2015, from Aubé-Michaud et al., 2017), and eight months with exclusive or notable snow input ; mean annual air temperature is -5. 2°C (1982-2012), with means of +9.5°C in July and -21.1°C in January. Annual precipitation is approximately 460 mm, with 45 % falling as snow. In Kangiqsualujjuaq, the SILA station of the Centre d'études Nordiques (CEN) indicates that winds come mostly from the north-north-west (Figure 2B) at an average speed of 19 km/h (during winter months, recording snowfall, the wind orientation varies from NNW to NNE). In general, the average daily wind speed is less than 49 km/h (99.8 % of the time), but hourly peaks of more than 50 km/h have been observed during the winter months [Aubé Michaud et al., 2017].

4

In the analysis of snow-avalanche hazards and associated risks, the knowledge of past events is highly depending on the occupation period and clear account of events (dates, location of the path, type of snow, runout distance, etc) [START_REF] Mcclung | The Avalanche Handbook[END_REF]Stethem et al, 2003], enabling differentiating the situation of risk from punctual catastrophic events [Dauphiné and Provitolo, 2003] ; in Kangiqsualujjuaq valley, human occupation has been only temporary before the settlement, the precise camp sites being not precisely known ; the village further recorded a constant demographic growth since its establishment. Wind has a strong role in the gathering of the snow that will be forming the snow avalanche : in Kangiqsualujjuaq, none of the known snowavalanche events originate from the usual NNW direction, but from the WSW to SW unusual blizzards (Figure 2B). A recent village with a growing population A historical perspective of Kangiqsualujjuaq 5 Kangiqsualujjuaq, which means "very large Bay", is the easternmost village in Nunavik and located in Akilasakallak Cove on the George River (Figure 1). Kangiqsualujjuaq was funded on a land that has been previously occupied by Pre-Dorset, Dorset and then by Thule peoples, as indicated by the archeological remains (Avataq Cultural Institute 2011). Kangiqsualujjuaq began as a series of temporary trading establishments along the George River estuary and the shore of the Ungava bay (Archives of Manitoba). The closest outspost was the Hudson's Bay Company outpost established in 1838 in the Illutaliviniq region. This outpost was named Fort Severight, and was located about 15 km south of Kangiqsualujjuaq. The Hudson Bay Company operated there from the first part of the 19 th century to the mid 20 th century as a salmon and seal fishery. Over 300 km north, along the shore of the Ungava Bay, the long lasting settlements of Port Burwell/Kiliniq has been occupied as a trading post or other posts from the end of the 19 th century to late 1970s.

6

The present-day location of Kangiqsualujjuaq was established in 1959 (Figure 3), when local Inuit created one of the first co-operatives in Northern Quebec to trade Arctic char with southern Canada (Archives of Manitoba). The site and its surroundings have long been occupied by temporary camps, as it is remembered by Tivi Etok from when he was a young boy [Weetaluktuk and Bryant, 2008]. However it is unknown whether these were summer or winter camps; these camps never rose to the level of a "village",

The development of Kangiqsualujjuaq and the threat of snow avalanches in a pe... while Killiniq was clearly identified as a large village well located on the migration route of marine mammals known by Inuit ancestors, "the oldest continuously inhabited site" [Tivi Etok in Weetaluktuk and Bryant, 2008] (Figure 1c). No evidence of continuous occupation of Kangiqsualujjuaq territory is found prior to 1959. From 1959 on, constructions were uninterrupted and included several public services (schools, government buildings, a health dispensary, a water supply system, a church, a municipal garage, a dump, gravel roads, a runway, sport infrastructures, etc) and numerous houses. Nowadays, Kangiqsualujjuaq is the starting point for visitors to the northeasternmost Canadian National Parks (Torngats Mountains and Kuururjuaq) and home of the park headquarters. A research station has also recently been established in the village ( 2018), managed by the University Laval and the community of Kangiqsualujjuaq. Regional Administration area. In Kangiqsualujjuaq, the population evolution has not been steady, although it has been always positive: the village never lost inhabitants, and its population grew from +3.5 to 58% from 1966-2016 (Figure 4B). The greatest increase in population occurred after the village gained its official status in 1980. Therefore, the projections indicate +6.91%, +9.13% and +11.82% growth based on regional baselines, the 1996-2015 average or the 2006-2015 average respectively [in [START_REF] Gerardin | Background report -Kangiqsualujjuaq master plan and zoning bylaw[END_REF].

8 Such a demographic boom in a Nordic village has led to the expansion of the territory used for housing, service buildings and transport infrastructures.

A spatial perspective of Kangiqsualujjuaq 9

As mentioned earlier, permanent settlement in Kangiqsualujjuaq, at its location by Akilasakallak cove, on the eastern bank of George River estuary, started with the establishment of a Co-op, in 1959. Only two families were settled there at that time [Weetaluktuk and Bryant, 2008], and peoples from small camps nearby moved there to form the beginning of a community. At that time, the village was composed of tents and only one log cabin. The first permanent buildings, such as the saw mill and the school, where erected on the lower land; prefabricated houses were erected from 1962 on the first terraces, as well as the church, as shown on Figure 5A. This location is sheltered from the prevailing winds coming from the southeast to the northwest (Figure 2B; Germain, 2016). In 1972 (Figure 5B), the village became more organized, with a clearer separation between houses and utility buildings. The newer development at that time was closer to the shore. The village grew below the mountain slope on the right bank of the river that drains the valley with more houses, stores, a Co-op and municipal office (1984, Figure 5C). From that point on, development occurred at the bottom of the valley (1991, Figure 5D) and further down the valley, with the extension of the school that had been previously relocated by the slope previously, and the construction of the new airport at its present location. The village then extended further within the valley, towards the north-west, and occupies the lower terraces north of Akilasakallak cove (2002, Figure 5E). This is the location for the new municipal garage and the drinking water supply. Figure 5F-2015 shows the further development on the western part of the village along the south-east facing slope, including new housing, the police building and the day-care building. The road system has significantly developed, including long portions of paved streets within the village and mixed gravel and paved roads towards the airport. There are also offset facilities such as the used-water basins and the dump.

10 From the beginning of the settlement in the 1960s to the present days, the spatial extension of the village has been significant, multiplying the number of buildings (residential and services) from few scattered buildings to become a village, with numerous streets and cross-roads. The original valley snout has been extensively used to develop the village over a greater area, by the northwestern bank of the cove. Later, the housing programs spread into the remaining areas in the valley, further from the cove. The impact of the 1999 snow avalanche on the spatial distribution of the village 12 The catastrophic avalanche hit the gymnasium in Kangiqsualujjuaq in the New Year's Eve 1999 at 01:45 [Stethem et al., 2003;[START_REF]Snow avalanche hazard assessment and risk management in northern Quebec, Canada[END_REF]. To celebrate the New Year, over 100 of the people were gathered in the gymnasium, which was located in the westernmost part of the school [START_REF] Campbell | Current and future snow avalanche threats and mitigation measures in Canada[END_REF]. At this time, the school and its gymnasium were the largest building close to the north-east facing slope (Figure 5D). The distance between the slope and the back of the gymnasium building was short, not exceeding 40 m [Stethem et al., 2003]. The snow passed through the gymnasium backwall by an additional 35 m, reaching a runout of approximately 75 m and extended along both sides of the gymnasium. It attained a width of 130 m in the deposited area. Due to the convex slope at the top of the rockwall (Figure 6), the slight leeward depression facilitated the gathering of blown snow in the windy environment, forming snow cornices [START_REF]The vulnerability of northern cities to weather-related hazards: case studies from the Province of Quebec, Eastern Canada[END_REF]. A 220 m wide snow avalanche, about 4 m thick, detached from a height of about 110 m a.s.l. on a slope 38° steep [Lied and Doomas, 2000;Stethem et al., 2003]. The southwest walls of the gymnasium (on which the main door was located), was breached with an impact pressure of about 10 kPa [START_REF] Schaerer | Avalanche hazard mitigation for Kangiqsualujjuaq, Québec[END_REF]. A blizzard had been in the area for several hours with western winds from 25-40 km/h following several days of blowing snow [START_REF] Schaerer | Avalanche hazard mitigation for Kangiqsualujjuaq, Québec[END_REF]; such winds were uncommon (Figure 2). The death toll was 9 dead (including 5 children)

The development of Kangiqsualujjuaq and the threat of snow avalanches in a pe...

Espace populations sociétés, 2020/3-2021/1 | 2021 with 25 people injured (about 1.5% and 4% of the village population at that time respectively). In black, a schematic reconstruction of the location of former buildings, which today corresponds to the exclusion zone, and known damaged buildings by snow-avalanche events prior to 1999; the white dashed lines represent the departure zones of those snow avalanches.

Photo: A. Decaulne.

13 Previously, other snow-avalanche events had been observed (Figure 6A), but these did not cause such extensive damage and harm to people. The dates for these events are not precise [START_REF] Schaerer | Avalanche hazard mitigation for Kangiqsualujjuaq, Québec[END_REF]Hétu, 2001]: the back wall of the school gymnasium was hit in March 1993, and the steps leading to the door were buried under 1.5 m of snow; two people were buried but extracted alive, although one was injured; the runout distance was 40 m.

14 None of these events had been formally reported, as the damage they caused was limited. After the 1993 event, however, the back of the school was considered unsafe by some of the inhabitants and an expert [Hétu, 2001] recommended that the snow cover be stabilized at the source-area and that the backdoor of the gymnasium be sealed.

15 Snow-avalanche deposits were observed at the foot of the slope in February 2000, during experts visit [Lied and Doomas, 2000], local people saying that such snow avalanches, however qualified as "not important" as provoking no damage, do happen yearly.

16 The devastating avalanche event in 1999 prompted the investigation of snow avalanches, highlighted the previous events that had not been formally reported, and prompted the Quebec Ministry of Public Safety to consult with various experts. It was recommended that the exclusion zone be further extended from 75 to 100 m [START_REF] Guilbault | Kangiqsualujjuaq, mitigating avalanche risk through land-use planning[END_REF]. This decision required the relocation of over 40 buildings, which mainly included houses but also a few business and services (a school, Co-op, church, etc.). The affected area also represented about 1/5 of the surface of the village at that time, requiring the Challenges ahead: conflicts between human needs and natural constrains 17 As the population of Kangiqsualujjuaq continues to grow, although at a slower pace that it did during the end of the 20th century, the need for housing also grows. The core of the village is already densely built. The main issue is to determine the location for new housing developments, along with associated infrastructure, such as transportation. The topography of the valley certainly represents a constraint, as slopes over 35° are prone to release slab snow avalanches of all sizes [START_REF] Mcclung | The Avalanche Handbook[END_REF]. Under specific snow conditions, when large amounts of windblown snow are falling over ice crusts formed into the previous snowpack, or during sudden thaw, slopes from 10 to 25° are prone to releasing infrequent slab avalanches, or loose wet snow avalanches. Slush flows can also be released on slopes of 5 to 25° [START_REF] Mcclung | The Avalanche Handbook[END_REF]. A wet snow avalanche will decelerate on slopes with a gradient lower than 10°. Rockfall cannot be excluded either, as observed during field work (Figure 7). Regarding the possible return period of the largest potential snow-avalanche event, the gathering of accurate data over 200 years at a specific location would be required to be 90% confident that a 100-year return period event had been documented within the chronology [START_REF] Mcclung | The Avalanche Handbook[END_REF]. This 200-years snow-avalanche history is not yet recorded in Kangiqsualujjuaq. Therefore, all areas along the southeastern slope are potentially at risk, and the greater exclusion zone established after the 1999 fatal snow avalanche is fully justified; however, the road to the airport lies very close to the length of the slope west of the village. Furthermore, the new housing development located along the gentler northwestern slope could also be affected by wet and slush snow avalanches as the slope gradients there are rarely below 10°.

18 The village is also located within a zone of abundant discontinuous permafrost, at the edge of the continuous permafrost zone [Allard et al., 2012]. The permafrost distribution and thickness, as well as the thickness of the active layer, vary considerably depending on the grain size and type of superficial sediments, snow-cover thickness, hydrological conditions and vegetation. In addition some zones of the valley are poorly drained (Figure 7).

19 Globally, the territory of the municipality of Kangiqsualujjuaq can be divided into three sectors according to its superficial sediment/substrate and topography (Figure 6) and associated reaction to permafrost thaw: bedrock outcrops, which can be divided into summit areas and slopes. These areas are exposed to dominant winds or the slope gradient prevents human occupation and infrastructure setting such as transportation corridors. The active layer is 4.5 to 6 m thick;

joints and cracks may contain ice. The rockmass may be unstable superficially where the slope gradient is steep [Aubé-Michaud et al., 2017].

marine deposits, composed of silt to fine sand, occupy the bottom of the valley. The active layer reaches 0.5 to 1.5 m; lots of segregation ice is present. The sediment subsides easily. In the center of the valley, between the airport and the village, there are abundant ice-rich poorly drained silt deposits. These thaw-sensitive deposits are unsuitable to construction, which limits the expansion of buildings in this area. These deposits of sandy fine sand are sensitive to thaw, creating wetland that are fed by surface water during frost-free seasons.

Unstable terrains during the thawing period are evidenced by the thermokarstic subsidence that affects and deforms the road. Furthermore, as in many places in the valley, several seasonal ice-core mounds were observed along the airport access road in a depressed area, in a fine sand deposit. Power poles along the road have been damaged because of this [L'Hérault et al., 2014;Aubé-Michaud et al., 2017].

glacial till, over 2 m thick in the lower parts, and less than 2 m thick at the apex, on the neighboring gentle slopes, are mostly present on the northern slope of the valley. The active layer is 2.5 to 3 m thick. This sediment is creeping and locally subsiding, i.e., unstable during the thawing period [Aubé-Michaud et al., 2017].

20 In such topographic and permafrost conditions, very little terrain is suitable for construction (Figure 7). The bedrock is made of granitic gneiss and would support buildings where the topography is not too uneven; there, the active layer has little effect on the surface stability. Sand and gravel deposits about 2 m thick over the bedrock could also be suitable, as the ice volume within such material remains below 10 m. Gravel and sand deposits over 2 m thick are also well drained, inhibiting important ice content, even with an active layer 2.5-3 m thick. Therefore, the most extensive areas suitable for further construction that could possibly support housing developments are (i) the site beside Akilasakallak cove where houses were transferred from the snow-avalanche exclusion zone after the 1999 avalanche, and (ii) the site beyond the airport. The first area mentioned still has some space left for new construction; the second one is located far from the village core, but already benefits from transportation infrastructures since part of it is used as a sand and gravel quarry. Other surfaces are not suitable for sustained development, requiring post-construction building adjustments and infrastructure maintenance. With reference to Figure 8, we observe that the whole valley bottom is excluded from the prone building surface, as marine deposits are unstable during the thawing period. From these profiles, one can observe the distribution of the three main surfaces within the valley, i.e., bedrock, glacial till and fine marine deposits. The indications regarding slope gradients also provide information on the potential slope processes, such as snow avalanches, rockfall and creep.

Elevation data extracted from LiDAR 2010, ministère de l'Énergie et des Ressources naturelles. The topography and permafrost situation in Kangiqsualujjuaq has forced the growing community to consider adaptive measures to reduce the risk posed to the population and to local infrastructures. A large part of the lowland area, at the bottom of the valley and along the coastline includes terrains that are unsuitable for construction. However, this is where most of the village is located and where the recent development has occurred. Snow-avalanche events forced the relocation of part of the village, and permafrost thawing associated with the deepening of the active layer required significant engineering design to sustain drainage and building stability. Areas suitable for construction are seldom found at low altitude and are mostly located west of the pass where the airport is located. There is greater exposure to wind here and the distance to the center of the village (and to facilities such as the school, medical center, stores, and municipal services) is much longer. Northern stakeholders have important choices to make, in collaboration with the residents (local observers) and scientists, about how to manage slope processes and respond to surface surveys that identify permafrost conditions at the local scale.

ABSTRACTS

Slope processes are active in the rolling plateau landscapes of Nunavik, Northern Québec, Canada. There are a few short but very steep escarpments in this region. On January 1 st , 1999 a powerful snow avalanche struck Kangiqsualujjuaq, one of the 14 Inuit villages in Nunavik. Nine people died and 25 were injured. This village and its surrounding are located within a glacial valley, in a periglacial environment. There is significant population growth, as well as in the other villages in Nunavik. As early as 1999-2000, there was a significant spatial reorganization of the village's infrastructures to avoid the impact of other snow-avalanche events. The main objective of this paper is to examine the village expansion in response to snow-avalanche process and population growth, within an area constrained with permafrost thawing and steep slopes.

From naturalist geomorphologic methods, written sources such as archive documents and aerial photographs, the results show that slopes above Kangiqsualujjuaq are prone to release snow avalanches during blizzards from uncommon directions, and that the newly built housing may be at risk in some places, due to the conjunction of snow avalanches and permafrost thawing.

Dans les paysages de plateaux ondulés du Nunavik, au nord du Québec, au Canada, les processus de versant sont actifs. Les escarpements sont rares et courts, mais sont significativement pentus.

Le 1 er janvier 1999, une avalanche mortelle a touché le village de Kangiqsualujjuaq, l'un des 14

Figure 1 -

 1 Figure 1 -Location map of Nunavik (A) in Northern Quebec (B) and Kangiqsualujjuaq surroundings (C), with National Parks in grey shade; a view of the village taken from the east-north-east at low tide (D).

  Photo A.[START_REF] Decaulne | Geomorphic evidence of slope dynamics on the Canadian Shield throughout the Holocene -a study from Lac à l'Eau-Claire, western Nunavik[END_REF]

  2

  Figure 2 -Climate data. A: temperature and precipitation in Kuujjuaq, from 2006-2015; B: origin and mean speed of daily winds in Kangiqsualujjuaq CEN station from 2005-2008 and winter months details (modified from Aubé Michaud et al., 2017).

Figure 3 -

 3 Figure3-Key dates in the history of Kangiqsualujjuaq (modified from[START_REF] Gerardin | Background report -Kangiqsualujjuaq master plan and zoning bylaw[END_REF] MAMOT 2015)

Figure 4 -

 4 Figure 4 -Demographics in Kangiqsualujjuaq, with raw data and 2036 projections (A) and population growth (B).

11

  With the opening of three National Parks in the surroundings of Kangiqsualujjuaq over the last 10 to 15 years (Torngats Mountains to the North, Kuururjuaq to the west and Ulittaniujalik by Pyramid Peak 115 km to the south), the village acquired the National Parks office and became the official gateway. Kangiqsualujjuaq is therefore welcoming many tourists to the area, and requires additional infrastructures (transportation, lodging, activities).

Figure 5 -

 5 Figure 5 -The spatial expansion of Kangiqsualujjuaq from 1967 to 2015, mapped from available aerial photographs.

Figure 6 -

 6 Figure 6 -Kangiqsualujjuaq in 2019.

  the "small store" was hit in March or April 1981 or 1982, catching one person who was not injured; the runout of the snow-avalanche from the foot of the slope was approximately 20 m; the back wall of the Co-op store was hit twice, in March orApril 1982April , 1984April or 1986: canoes : canoes were damaged in one of the two snow avalanches, and one person was caught in the other one, without sustaining an injury; the runout distance to the Co-op was about 20 m; the back wall of the Youth Centre was hit in March or April 1984 or 1986; no one was injured, although the runout distance to the building is the shortest, at only 14 m;

  Kangiqsualujjuaq and the threat of snow avalanches in a pe... Espace populations sociétés, 2020/3-2021/1 | 2021 designation of new building locations along Akilasakallak cove and the northern part of the valley (Figures 5 E-G and 6).

••

  The development of Kangiqsualujjuaq and the threat of snow avalanches in a pe... Espace populations sociétés, 2020/3-2021/1 | 2021

•Figure 7 -

 7 Figure 7 -Cross profiles of Kangiqsualujjuaq valley.

Figure 8 -

 8 Figure 8 -Location of the suitable construction areas on surfaces where permafrost thawing and active layer fluctuations have little effect on the stability of the surface, away from active slope processes areas (modified from Aubé-Michaud et al., 2017).

  photographies aériennes diachroniques, les résultats montrent d'une part que les versants qui dominent Kangiqsualujjuaq sont favorables au déclenchement d'avalanches durant les blizzards de directions inhabituelles, et d'autre part que les logements récemment construits peuvent dans certains cas se trouver également en situation de risque, en raison de la conjonction des avalanches et de la fonte du pergélisol.

  The development of Kangiqsualujjuaq and the threat of snow avalanches in a pe... STETHEM C., JAMIESON B., SCHAERER P., LIVERMAN D., GERMAIN D., WALKER S., 2003, Snow avalanche hazard in Canada -a review. Natural Hazards, 28, pp. 487-515. VEILLEUX S., BHIRY N., DECAULNE A., 2020, Talus slope characterization in Tasiapik Valley (subarctic Québec): evidence of past and present slope processes. Geomorphology Doi : 10.1016/ j.geomorph.2019.106911

	nordiques, Université Laval, 68 p.
	ARCHIVES OF MANITOBA, http://pam.minisisinc.com/scripts/mwimain.dll/144/
	PAM_AUTHORITY/AUTH_DESC_DET_REP/SISN%20115?sessionsearch
	AVATAQ CULTURAL INSTITUTE, 2011, http://www.avataq.qc.ca/en/Institute/Departments/
	Archaeology/Discovering-Archaeology/Arctic-Chronology.
	BÉRUBÉ J., 2000, Rapport d'enquête publique suite à l'avalanche de 1999 à Kangiqsualujjuaq. Rapport du
	bureau du coroner, 66 p.
	BHIRY N., DECAULNE A., BOURGON-DESROCHES M., 2019, Development of a subarctic peatland
	linked to slope dynamics at Wiyasakami Lake (Nunavik, Canada). The Holocene, 29, 9,
	pp. 1459-1467.

BIBLIOGRAPHY ALLARD M., LEMAY M., BARRETTE C., L'HÉRAULT E., SARRAZIN D., 2012, Permafrost and climate change in Nunavik and Nunatsiavut: Importance for municipal and transportation infrastructures. In Allard M. and Lemay M. (Eds): Nunavik and Nunatsiavut: From science to policy. An Integrated Regional Impact Study (IRIS) of climate change and modernization. ArcticNet Inc., Québec, Canada, pp. 171-197. AUBE-MICHAUD S., ALLARD M., L'HERAULTE., MATHON-DUFOUR V., 2017, Identification des risques actuels et appréhendés sur le territoire des communautés du Nunavik en fonction des changements climatiques -phase 1 : Document de synthèse, communauté de Kangiqsualujjuaq. Rapport final réalisé pour le compte du Ministère de la Sécurité publique (Gouvernement du Québec), Centre d'études Espace populations sociétés, 2020/3-2021/1 | 2021

The development of Kangiqsualujjuaq and the threat of snow avalanches in a pe... Espace populations sociétés, 2020/3-2021/1 | 2021

This work is supported by the French National Research Agency through the Labex DRIIHM program « Investissements d'avenir » referenced ANR-11-LABX-0010. It is also supported by the French Polar institute IPEV (funding #1148 DeSiGN), and the Natural Sciences and Engineering Research Council of Canada (NSERC). We also thank the council of Kangiqsualujjuaq and the inhabitants of the village, our partners at the Kativik Regional Government and the Makivik corporation for their valuable collaboration and the Centre d'études nordiques for logistical support. We also thank Emmanuel L'Hérault for providing elevation data to build up the longitudinal profiles. Kangiqsualujjuaq and the threat of snow avalanches in a pe...

AUTHORS


Marco Aggravi
email: marco.aggravi@irisa.fr

L Estima

Alexandre Krupa
email: alexandre.krupa@inria.fr

Sarthak Misra
email: s.misra@utwente.nl

Claudio Pacchierotti
email: claudio.pacchierotti@irisa.fr

Daniel A L Estima

Haptic

Keywords:

HAL is

I. INTRODUCTION

M OST robotic systems for flexible needle insertion still lack effective feedback capabilities [START_REF] Abayazid | Experimental evaluation of co-manipulated ultrasoundguided flexible needle steering[END_REF]- [START_REF] Pacchierotti | Teleoperation of steerable flexible needles by combining kinesthetic and vibratory feedback[END_REF]. In this respect, real-time visualization of the needle pose has been proven effective in increasing the clinician's performance in many of these scenarios [START_REF] Gerovichev | The effect of visual and haptic feedback on manual and teleoperated needle insertion[END_REF], using, e.g., magnetic resonance imaging (MRI) [START_REF] Moreira | Tele-operated mri-guided needle insertion for prostate interventions[END_REF], computed tomography (CT) [START_REF] Ben-David | Evaluation of a ct-guided robotic system for precise percutaneous needle insertion[END_REF], or ultrasound (US) [START_REF] Abayazid | Experimental evaluation of ultrasound-guided 3d needle steering in biological tissue[END_REF]. Most teleoperated robotic systems for needle insertion use one of these imaging modalities to provide visual feedback on the needle during the insertion. For example, Chevrie et al. implemented semi-automatic teleoperation of a flexible needle based on US tracking [START_REF] Chevrie | Flexible needle steering using ultrasound visual servoing[END_REF], [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF]. It uses an iterative algorithm to fit a polynomial curve to the needle shape, while a 3-dimensional (3D) version of a Start Algorithm [START_REF] Friedland | Automatic ventricular cavity boundary detection from sequential ultrasound images using simulated annealing[END_REF] tracks a spherical target. In [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF], the operator controls the needle tip velocity using a grounded haptic interface, while being provided with visual feedback via three orthogonal US slices that intersect with the needle tip and kinesthetic navigation feedback.

Haptic feedback is another important piece of feedback in robotic teleoperation, which has been proven effective for providing both navigation [START_REF] Abayazid | Experimental evaluation of co-manipulated ultrasoundguided flexible needle steering[END_REF], [START_REF] Pacchierotti | Teleoperation of steerable flexible needles by combining kinesthetic and vibratory feedback[END_REF], [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF], [START_REF] Selvaggio | Haptic-based shared-control methods for a dual-arm system[END_REF], [START_REF] Abi-Farraj | A haptic shared-control architecture for guided multi-target robotic grasping[END_REF] and environmental information [START_REF] Meli | Experimental evaluation of magnified haptic feedback for robot-assisted needle insertion and palpation[END_REF], [START_REF] Chinello | A three revolute-revolute-spherical wearable fingertip cutaneous device for stiffness rendering[END_REF]- [START_REF] Pacchierotti | Cutaneous device for teleoperated needle insertion[END_REF]. For example, Pacchierotti et al. presented a haptic teleoperation system for flexible needles. It enables clinicians to maneuver the surgical tool while providing them with navigation guidance through kinesthetic and vibratory force feedback [START_REF] Pacchierotti | Teleoperation of steerable flexible needles by combining kinesthetic and vibratory feedback[END_REF]. Majewicz and Okamura used a haptic interface to command the desired position of the needle tip in Cartesian space and provide force feedback representing kinematic constraints and the position error of the robotic system [START_REF] Majewicz | Cartesian and joint space teleoperation for nonholonomic steerable needles[END_REF]. Seifabadi et al. developed a 5-degreesof-freedom (5-DoF) parallel pneumatically actuated modular robot for teleoperated prostate biopsy under MRI [START_REF] Seifabadi | Robotic system for mri-guided prostate biopsy: feasibility of teleoperated needle insertion and ex vivo phantom study[END_REF].

However, it is rare to find teleoperation systems conveying feedback about the needle tip cutting/penetration force, which is often considered to be important [START_REF] Meli | Experimental evaluation of magnified haptic feedback for robot-assisted needle insertion and palpation[END_REF]. This limitation is mostly due to the difficulty in isolating this force from other undesired components (e.g., friction). In this respect, De Lorenzo et al. proposed a coaxial needle insertion assistant that isolates the cutting force and provides it to the operator with kinesthetic feedback. A smaller inner needle is inserted inside a larger outer one. The outer needle covers the inner one, so as to prevent any friction force from acting on the latter. Then, the tip of the inner needle is left uncovered, so that it can penetrate the tissue while sensing the cutting force only [START_REF] Lorenzo | Experimental evaluation of a coaxial needle insertion assistant with enhanced force feedback[END_REF]. Similarly, Elayaperumal et al. [START_REF] Elayaperumal | Mr-compatible biopsy needle with enhanced tip force sensing[END_REF] and Khan et al. [START_REF] Khan | Force sensing in continuum manipulators using fiber bragg grating sensors[END_REF] used Fiber Bragg Grating (FBG) sensors to enhance force and shape sensitivity of medical instruments. Other approaches use external sensors, combined with needle and soft tissue models [START_REF] Okamura | Force modeling for needle insertion into soft tissue[END_REF], [START_REF] Van De Berg | The influence of tip shape on bending force during needle insertion[END_REF].

This paper proposes a haptic-enabled robot-assisted teleoperation system for the insertion of flexible needles in soft tissue. It enables intuitive steering of the surgical tool while providing rich feedback information on needle-tissue interaction and navigation toward the target objective. Fig. 1 shows the system, while Fig. 2 describes how each relevant part is interconnected.

With respect to other works on the topic, the proposed system provides an unprecedented amount of environmental and guidance information, conveyed combining grounded (kinesthetic), ungrounded (vibrotactile), and visual feedback. For example, Chevrie et al. also presented an ultrasoundguided teleoperation system for flexible needles, but it only provides navigation feedback through kinesthetic stimuli [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF]. The contributions of the paper can be summarized as follows:

• design of a real-time 3D friction force estimation algorithm, so as to provide feedback about the force applied by the needle tip when cutting/penetrating the tissue; • design of rich vibrotactile-kinesthetic haptic feedback, so as to simultaneously convey navigation and cutting force information in a complete yet intuitive manner; • integrate the proposed methods with a state-of-the-art ultrasound-driven needle steering algorithm, to attain an effective robot-assisted needle insertion system; • carry out an extensive human subjects evaluation, aimed at finding the best feedback rendering technique as well as evaluate the system's overall performance (precision, repeatability, completion time, user's preference).

II. MATERIALS AND METHODS

A. Experimental setup

The remote system is shown in Fig. 1a. It consists of a Chiba biopsy needle (Angiotech MCN2208, CA) mounted on a 6-DoF force/torque sensor (ATI Nano 43, USA), which is in turn attached to a 6-DoF robotic manipulator (Omron Viper S650 Adept, JP). The needle length is 126 mm, with internal and external diameters of 0.48 mm and 0.70 mm, respectively. The needle's tracking is guaranteed by a 3D ultrasound (US) scanner (Analogic SonixTOUCH, USA) with a 4DC7 3/40 motorized 3D US probe. The 3D volumes are reconstructed in Cartesian coordinates after a 3 × 3 × 3 median filter is applied to remove some noise in the pre-scan data [START_REF] Lee | Intensity-based visual servoing for non-rigid motion compensation of soft tissue structures due to physiological motion using 4d ultrasound[END_REF]. The steering control and the US-guided tracking algorithms are described in a previous work of our group [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF] (see Sec. II-B). A gelatine phantom is prepared following the procedure in [START_REF] Van De Berg | The influence of tip shape on bending force during needle insertion[END_REF].

The local system is shown in Fig. 1b. It consists of a grounded haptic interface (Force Dimension Omega.6, CH), a custom wearable vibrating armband, and a screen. The Omega.6 is used as an impedance haptic interface: we measure the position of its end-effector/handle, controlled by the human operator, to steer the needle; at the same time, through the same handle, we provide the operator with kinesthetic feedback from the remote environment. The vibrating armband comprises four actuators, positioned evenly around the arm at 90 degrees from each other [START_REF] Aggravi | Design and evaluation of a wearable haptic device for skin stretch, pressure, and vibrotactile stimuli[END_REF]. Finally, the screen shows the current ultrasound images and a real-time 3D reconstruction of the needle structure with respect to the target. The overall architecture of our teleoperated needle insertion system is illustrated in Fig. 2.

B. Needle steering algorithm

The needle steering algorithm translates the velocity commanded by the user for the needle tip, v h , into a 6-DoF velocity command for the robotic manipulator which houses the needle, v n . To do so, we use a Jacobian matrix J s obtained with the redundancy formalism of the task-function framework [START_REF] Slotine | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]. This matrix links the 6-DoF velocity of the needle tip Cartesian frame {T} to that of the needle base Cartesian frame {B} (see Figs. 2 and3). This control acts similarly to a velocitycontrol teleoperation (displacements at the local side with respect to the initial position are translated into velocities of the robotic end-effector), modified to consider task priorities as in the stack-of-tasks framework [START_REF] Slotine | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]. In this respect, we consider three tasks: i) ensure that the needle base stays aligned with the insertion point, limiting the lateral tissue deformation at the entry point; ii) have the direction of the tip velocity oriented toward the target (this task generates our navigation feedback, see Sec. III-A); and iii) ensure that the circular tip trajectory induced by the bevel is oriented in the desired direction, minimizing the needle bending.

In this work, we use the same framework presented in [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF], [START_REF] Chevrie | Flexible needle steering in moving biological tissue with motion compensation using ultrasound and force feedback[END_REF], [START_REF] Chevrie | Online prediction of needle shape deformation in moving soft tissues from visual feedback[END_REF] for tracking the needle shaft from 3D ultrasound and determining the Jacobian that links the needle tip velocity to the base velocity. Moreover, as discussed in Sec. III-A1 and shown in Fig. 2, this framework is also used to compute the guiding reference f g , which helps the human operator driving the needle toward the target.

C. Friction Force Estimation

As we are interested in rendering the insertion/cutting force at the needle tip but our sensor is only placed at its base, we need to estimate (and then filter out) the friction force acting along the needle shaft.

The total force applied on the needle f s (l) ∈ R 3 at its base is registered by our F/T sensor. After we compensate for the needle (and its support) gravity contributions, force f s (l) is due to the friction between the needle shaft and the tissue f f ∈ R 3 as well as to the cutting force of the needle tip f t ∈ R 3 [START_REF] Okamura | Force modeling for needle insertion into soft tissue[END_REF] (see Fig. 3),

f s (l) = f f (l) + f t (l), (1)
where l ∈ R + is the length of the needle inside the tissue. Due to the needle bending, the friction force along the shaft has components in three directions (3D). Roesthuis et al. computed these components in a 2D plane using the needle slope as the angle between the needle's centre axis and the coordinate

Control

Remote system Operator Fig. 2. Haptic teleoperation system. The operator uses a grounded kinesthetic device and a wearable vibrotactile armband (left). A needle steering algorithm (center) translates the human operator's velocity commands v h into a velocity vn for the 6-DoF Viper robot (Omron, JP), which houses a 6-axis force/torque (FT) sensor (ATI, USA) and the needle (right). During the insertion, a friction force estimation algorithm filters the FT sensed forces f s and estimates those due to the penetration of the needle tip into the tissue, f t . A 3D ultrasound (US) probe enables the tracking of the needle x tip and target xtarg at runtime, while a steering algorithm calculates the best needle tip direction to reach the target objective, d n tip . The user receives navigation guidance (G) and needle cutting force (C) using a combination of kinesthetic (subscript K, forces fg and f cf ) and vibrotactile (subscript V, vibrating motors frequencies mg and m cf ) feedback. On the screen, the operator sees the ultrasound images and a real-time reconstruction of the needle pose. frame [START_REF] Roesthuis | Mechanics of needle-tissue interaction[END_REF]. Extending this model to 3D, we can estimate the friction force by using the tangent to the needle shaft

f f (l) = L L-l f f t(x b (c))dc , (2)
where x b (c) ∈ R 3 is the position of a point of the needle at the curvilinear coordinate c, t(x b (c)) ∈ R 3 corresponds to the normalized tangent vector at the location x b (c), f f ∈ R is defined as friction per unit inserted needle length [START_REF] Roesthuis | Mechanics of needle-tissue interaction[END_REF], l is again the inserted length, and L is the total needle length. By using (2), the needle 3D reconstruction, and the inserted length estimated from the US tracking [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF], we can obtain the total friction force f f (l) over the needle inside the tissue, l. However, f f depends on several factors, such as the tissue properties and insertion velocity, and, in viscoelastic tissues, it will increase as the relative velocity of the needle with respect to the tissue increases [START_REF] Roesthuis | Mechanics of needle-tissue interaction[END_REF]. To ease our estimation, we take advantage of the fact that in viscoelestic tissues the friction force has the same amplitude if the needle is being inserted or removed, as long as the needle velocity is the same [START_REF] Asadian | Dynamics of translational friction in needle-tissue interaction during needle insertion[END_REF]. For this reason, when removing the needle from the soft tissue, we can assume that f t (l) in ( 1) is equal to zero, thus leaving us with f s (l) = f f (l).

We approximate the non-straight configuration of the inserted needle to a series of smaller straight segments. Then, we can calculate each segment's tangent t(x b (c)) and corresponding f f values. Summing the friction contribution of all segments provides us with the final friction force applied on the needle shaft, f s (l) [START_REF] Roesthuis | Mechanics of needle-tissue interaction[END_REF]. Since the curvature of the needle affects friction [START_REF] Roesthuis | Mechanics of needle-tissue interaction[END_REF], a correct estimation of the friction force also relies on the length of the considered segments (more segments equals a better estimation but also a longer computation time). In our case, we used ten segments that start really small and grow as the needle penetrates the gelatine, reaching a maximum length of 13 mm when the needle is fully inserted. In this respect, it would be useful to adjust the number of considered segments to the uncertainty of the current estimation, so as to always guarantee a minimum estimation quality.

To estimate f f , we collected frictional data at different insertion velocities, and then we used such recording to Fig. 3. The robotic insertion system. The needle is inserted into the soft tissue phantom at velocity vn. The force/torque (FT) sensor registers the total force applied to the needle fs, which is composed of the cutting force at the needle tip ft and the friction force along its shaft f f . The ultrasound (US) tracks the pose of the needle tip x tip and the target xtarg at run-time. The steering algorithm calculates the projection of xtarg onto the plane perpendicular to the needle P ⊥ , i.e., x P targ , so as to provide navigation information minimizing the distance d n tip between x tip and x P targ . Needle insertion guidance toward the target, cutting force, the US images, and a reconstruction on the needle are provided to the user through a combination of visual, kinesthetic, and vibrotactile feedback. perform a linear regression for estimating the friction force per unit length f f for any insertion velocity. Of course, this approach can be easily extended to consider more complex friction models [START_REF] Armstrong-Hélouvry | A survey of models, analysis tools and compensation methods for the control of machines with friction[END_REF]. As Chevrie et al. [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF] considered 0.005 m/s as the maximum insertion velocity, we included these five velocities in our analysis:

v = [v 1 , . . . , v 5 ]
T = [0.001, 0.002, 0.003, 0.004, 0.005] T m/s. For each velocity v i , i = 1, . . . , 5, (that corresponds to the velocity of the needle along the z-axis of B, i.e., v z n ), we recorded the F/T sensor readings for three 50-mm-long insertions inside the soft-tissue phantom. Then, for each trial, we performed a linear regression to obtain the slope coefficient, and we averaged the results for each velocity. We only considered the data registered during the retraction of the needle. This approach provides us with a set of five f f estimates f e f = [ ff,v1 , . . . , ff,v5 ] T , one for each considered velocity. We can now estimate at runtime the current f f , depending only on the current needle axial velocity v z n ,

f f (v z n ) =    v z n ff,v 1 v1 if v z n < v 1 , ff,vi + ff,v i+1 -ff,v i vi+1-vi (v z n -v i ) if v i < v z n < v i+1
.

(3) The upper term represents a linear regression between v 1 and zero velocity, while the lower term applies a linear regression between the two closest velocity estimates. Figure 4a shows an example of insertion with constant velocity along the needle z-axis (needle shaft) v z n of ±0.005 m/s. It highlights the registered readings along the z-axis of the F/T sensor f z s (left ordinate axis, light blue) and the needle velocity along the direction of its shaft v z n (right ordinate axis, orange), over time. It is visible how f z s linearly decreases with the retraction (positive sign) or insertion (negative sign) of the needle. Figure 4b shows a detail (red box in Fig. 4a) of the force during the retraction of the needle, for all the considered velocities. Fig. 4c reports all the obtained f f coefficients. As expected, faster velocities result in higher force variations.

These data are used to filter out the friction force acting on the needle shaft from the total force registered by the F/T sensor. Doing so, we can provide the user with an estimate of the insertion/cutting force at the needle tip. Of course, our feedback rendering algorithms described below can use of any similar friction estimation algorithm, e.g., [START_REF] Lorenzo | Experimental evaluation of a coaxial needle insertion assistant with enhanced force feedback[END_REF], [START_REF] Elayaperumal | Mr-compatible biopsy needle with enhanced tip force sensing[END_REF], [START_REF] Okamura | Force modeling for needle insertion into soft tissue[END_REF].

III. FEEDBACK INFORMATION

The user drives the needle inside the soft tissue using a grounded haptic interface, as described in Sec. II-B. The system provides him or her with two types of feedback information: needle direction guidance (Sec. III-A) and cutting force rendering (Sec. III-B). The former helps the user steering the needle toward the given objective (e.g., the area to treat), while the latter provides information about the estimated insertion/cutting force at the needle tip. Both pieces of information are important, and we evaluate how it is best to provide them: through kinesthetic feedback (Secs. III-A1 and III-B1), delivered by the grounded haptic interface, or cutaneous feedback (Secs. III-A2 and III-B2), delivered by a vibrotactile armband.

A. Needle direction guidance (G)

Let us define x targ ∈ R 3 as the position of the target point (see Fig. 3), x tip ∈ R 3 that of the needle tip, n dir ∈ R 3 as the normalized direction vector of the needle, and P ⊥ as the plane originating in x tip and with orthogonal vector n dir . Let us also define x P targ as the the orthogonal projection of x targ onto P ⊥ , and d tip = x P targx tip as the vector between the target point and the needle position. Our guiding problem, i.e., driving the user along the direction toward the target, becomes that of minimizing d tip . In other words, the closer the needle tip position x tip is to x P targ , the more aligned it is toward the target. In order to provide this guiding information with respect to the operator, as to make it easy and intuitive to understand, we transform the projected target position, where R n us is the rotation matrix bringing x P targ from being defined in W to T. Matrix R n us ∈ SO(3), indicating the rotation between the US frame and the needle tip, is available from the steering algorithm of [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF]. With this transformation, the three-dimensional distance d tip becomes two-dimensional with respect to the needle tip T, since x P,n targ lies in the same plane P ⊥ originated by x tip , i.e., d n tip = [ d n,x tip , d n,z tip ] T . Aligning the needle tip to point toward the target reduces d n tip to zero, since the origins of P ⊥ and x P,n targ coincide. The system can provide this guiding information using either 1) kinesthetic feedback or 2) vibrotactile cutaneous feedback.

x P,n targ = R n us x P targ , 80
1) Needle direction guidance using kinesthetic feedback (G K ): The feedback provided to the operator consists of a kinesthetic force f g ∈ R 3 applied by the Omega.6 haptic interface,

f g = -K g (x h -x h,0 ). (4)
where x h and x h,0 are the current and rest position of the Omega.6 handle, respectively, and K g ∈ R 3×3 is an anisotropic stiffness characterized by low stiffness values along the axis aligned with the needle's target direction, i.e., that minimizing d n tip , and high stiffness values along other axis [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF]. In this condition, users feel a low resistance when moving the haptic interface along the direction toward the target, while they feel a strong resistance when moving in any other direction. This is the only type of haptic feedback conveyed by Chevrie et al. [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF] during their teleoperation. Vibrotactile sensations guide the user toward minimizing d tip , i.e., the distance between the needle tip x tip and the projection of the target onto the plane perpendicular to the needle P ⊥ , i.e., x P targ . In this example, as t d < d n,x tip < tu, t d < d n,z tip < tu, -d n,x tip < t d , and -d n,z tip < t d , we calculate m(1), m(2) using the middle term in (5) and m(3), m(4) using the first term.

2) Needle direction guidance using vibrotactile feedback (G V ):

We map d n tip onto the four vibrating motors of our wearable haptic armband. As the motors are positioned evenly around the arm at 90 degrees from each other, we can easily provide navigation cues along the transverse plane cutting the arm.

To do so, let us introduce the step-like activation function m(x):

m(x) =        0 if x < t d x m + (x M -x m ) (x -t d ) (t u -t d ) if t d ≤ x ≤ t u x M if x > t u , (5)
where t d and t u indicate the length of the step. Function m(x) is zero until x reaches the lower limit t d , then linearly grows from x m to x M between t d and t u , and it finally plateaus at x M after that. Let us define m = [m 1 , m 2 , m 3 , m 4 ] as the commanded vibrating motors frequencies (see Fig. 1). The minimum and maximum perceived activation frequencies of our motors is 60 and 280 Hz, respectively, [START_REF] Scheggi | Cooperative navigation for mixed human-robot teams using haptic feedback[END_REF]. By considering this, and by setting t d = 0.002 m and t u = 0.010 m as activation thresholds, our mapping is described as follows

m g =     m 1 m 2 m 3 m 4     =     m(d n,x tip ) m(d n,z tip ) m(-d n,x tip ) m(-d n,z tip )     . (6)
Fig. 5 shows a representative example of this approach with d n tip = [ 3.8, 9.6 ] T mm. To reach the target, the operator needs to move the needle from x tip toward x P targ . The system conveys navigation information minimizing d n tip , i.e., since t d < d n,x tip < t u and t d < d n,z tip < t u , we use both times the middle term in [START_REF] Moreira | Tele-operated mri-guided needle insertion for prostate interventions[END_REF]. This leads to the activation of two motors in the armband, guiding the user to move the needle tip toward x P targ .

B. Cutting force rendering (C)

The second piece of information we want to provide is the cutting/insertion force at the needle tip, estimated with the help of the friction force estimation algorithm described in Sec. II-C. Again, the system can provide this information using either 1) kinesthetic feedback or 2) vibrotactile cutaneous feedback.

1) Cutting force rendering with kinesthetic feedback (C K ):

In Sec. III-A1, we use the anisotropic stiffness of matrix K cf to keep the user along the right direction of insertion. This matrix can be modified as to provide the operator with cutting force feedback instead of guiding directions. To do so, recalling (1), let us define ft = f t /f M ∈ R 3 as the adjusted cutting force at the needle tip, f M ∈ R + as the maximum force that the haptic device can rendered (12 N for the Omega.6), and

K h =   k h 0 0 0 k h 0 0 0 k h  
as a positive diagonal isotropic stiffness matrix (we consider k h = 300 N/m in our experiments). By setting

K cf = K h + diag(K h ft ), (7)
we are able to provide information about f t with the same approach used in Sec. III-A1, i.e., K cf is an anisotropic stiffness matrix with higher stiffness along the directions spanned by f t . Finally, as for Sec. III-A1, the feedback provided to the operator consists in a kinesthetic force f cf ∈ R 3 applied by the Omega.6 haptic interface,

f cf = -K cf (x h -x h,0 ). ( 8
)
2) Cutting force rendering with vibrotactile feedback (C V ): We use the same step-like activation function described in [START_REF] Moreira | Tele-operated mri-guided needle insertion for prostate interventions[END_REF]. However, this time we are interested in providing information about the force variation along the insertion axis f z t (see Fig. 1 and Fig. 3). By setting again x m = 60 Hz and x M = 280 Hz as the minimum and maximum activation frequencies of our motors, with t d = 0.25 N and t u = 0.75 N activation thresholds, our mapping is described as follows

m cf = m 1 , . . . , m 4 T = m(f z t ), . . . , m(f z t ) T . ( 9
)
In this condition, users feel no vibrations whenever the estimated cutting force stays below t d = 0.25 N. Then, all motors starts to vibrate with increasing intensity until the force reaches t u = 0.75 N. To avoid spurious vibrations due to noise in the estimation of the force, we used a moving average filter to smooth f z t before computing m(f z t ) (window size = 10).

IV. EXPERIMENTAL EVALUATION

To evaluate the effectiveness of our teleoperated system, we carried out a human subjects study. A video is available as supplemental material and at https://youtu.be/IXZ-n4hzLFM.

A. Experimental setup and task

We considered the robotic teleoperation system shown in Fig. 1 and described in Sec. II-A. The needle is inserted into a soft-tissue phantom made of gelatine mixture, to which silica powder is added to mimic the acoustic scattering of human tissue [START_REF] Abayazid | Experimental evaluation of co-manipulated ultrasoundguided flexible needle steering[END_REF]. A small elliptical target, i.e., a 3-cm-long olive, is embedded in the gelatine and it is placed at around 9 cm from the surface where the needle is inserted (see also Fig. 3).

Participants were asked to steer the needle inside the gelatine until they punctured the outer surface of the target. They were provided with a combination of navigation and needle tip force feedback according to the condition being considered (see Sec. III). Before the beginning of each insertion, the needle was perpendicularly pre-inserted into the gelatine for 8 mm, so as to be able to initialize the US tracking. The point of insertion of the needle was randomly chosen among seven options (see video). This approach enabled us to use one gelatine per participant and ensure that the needle had to be always steered along a non-straight path to reach the target.

B. Subjects

Thirteen participants took part in the experiment, including 2 women and 11 men (age 22-36 years old). One practice trial per feedback condition was allowed. 7 had previous experience with using haptic interfaces. Each participant tested one time each of the seven combinations of devices, for a total of 78 insertion trials. Users were asked to complete the task as precisely as possible taking into account the feedback received.

C. Feedback modalities

We considered seven different feedback modalities: G N C N neither needle direction guidance nor cutting force rendering is provided to the user; G N C K cutting force rendering only, provided by the Omega.6 grounded kinesthetic interface as detailed in Sec. III-B1; G K C N needle direction guidance only, provided by the Omega.6 kinesthetic interface as detailed in Sec. III-A1; G N C V cutting force rendering only, provided by the wearable vibrotactile armband as detailed in Sec. III-B2; G V C N needle direction guidance only, provided by the wearable vibrotactile armband as detailed in Sec. III-A2; G V C K needle direction guidance provided by the vibrotactile armband (Sec. III-A2) and cutting force rendering provided by the Omega.6 kinesthetic interface (Sec. III-B1); G K C V needle direction guidance provided by the Omega. [START_REF] Ben-David | Evaluation of a ct-guided robotic system for precise percutaneous needle insertion[END_REF] kinesthetic interface (Sec. III-A1) and cutting force rendering provided by the vibrotactile armband (Sec. III-B2); Visual feedback (US images and reconstruction of the needle) is always provided to the users through a screen placed in front of the grounded interface (see Fig. 6 and video). Fig. 7 shows a representative insertion in condition G K C V . The top figure shows the profile of the estimated cutting force f n t , while the bottom figure shows the evolution of the needle velocities during the insertion. When the needle punctures the olive, the force rises and activates the vibrotactile feedback (in blue).

D. Results

Figure 7 shows the cutting force profile for one representative insertion in condition G N C V . The blue-shaded area shows when the armband was activated to inform the operator about the cutting force. The targeting accuracy for this trial was 0.8 mm, which is higher than that recently obtained by Chevrie et al. [START_REF] Chevrie | Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance[END_REF] for a similar setup and task. However, Chevrie et al. only provided navigation information through kinesthetic feedback, similarly to what we do in G K . This promising result is confirmed by the analysis in the remainder of this Section.

As a measure of performance, we registered (i) the targeting accuracy (how close the needle tip is to the target), i.e., x tipx targ at the end of the task; (ii) the percentage of trials the needle reached the target; (iii) the insertion time, calculated as the time between the first movement of the needle and the end of the task; (iv) the average insertion speed, calculated as the integral of the end-effector velocities along the insertion direction; (v-vi) the perceived effectiveness as registered by the users, evaluated using subjective questionnaires.

To compare the different metrics, we ran one-way repeatedmeasures ANOVA tests (a= 0.05). All data passed the Shapiro-Wilk normality test. Since not all combinations between needle direction guidance (G N , G V , G K ) and cutting force rendering (C N , C V , C K ) are tested, and interactions between variables are expected, a two-way repeated measures design is not appropriate in this case. Figure 8a shows the (i) targeting accuracy. Mauchly's Test of Sphericity indicated that the assumption of sphericity had been violated (χ 2 (20) = 63.118, p < 0.001). The ANOVA test with a Greenhouse-Geisser correction revealed a statistically significant change in the task accuracy between conditions (F(1.819, 21.826) = 32.219, p < 0.001). Results of the post hoc analysis with Bonferroni adjustments are shown in Table I. As a complement, Figure 8b shows the (ii) percentage of successful insertions, i.e., when the needle punctured the target. Figure 8c shows the (iii) insertion time. The collected data passed the Mauchly's Test of Sphericity. The ANOVA test revealed no statistically significant change in the insertion time between conditions (F(6, 72) = 0.345, p > 0.05). Figure 8d shows the (iv) insertion speed. The collected data passed the Mauchly's Test of Sphericity. The one-way ANOVA test revealed no statistically significant change in the insertion speed between conditions (F(6, 72) = 0.764, p > 0.05). In addition to the objective evaluation reported above, we also measured users' experience. Immediately after the experiment, subjects were asked to report the effectiveness of (v) each one of the seven feedback condition, of (vi) kinesthetic and vibrotactile feedback in providing either needle direction guidance or cutting force. Answers were registered using scales from 0 (not effective at all) to 10 (very effective). Figure 9a shows the perceived effectiveness for the seven experimental conditions. Mauchly's Test of Sphericity indicated that the assumption of sphericity had been violated (χ 2 (20) = 50.997, p < 0.001). The one-way ANOVA test with a Greenhouse-Geisser correction revealed a statistically significant change in the perceived effectiveness between conditions (F(2.612, 31.339) = 36.013, p < 0.001). Results of the post hoc analysis with Bonferroni adjustments are shown in Table I. Figure 9b shows the perceived effectiveness for kinesthetic and vibrotactile feedback in providing either needle direction guidance or cutting force (C V vs. C K and G V vs. G K ). Two paired samples t-test determined that this metric differed statistically significantly between conditions (C V vs. C K : t(12) = 6.052, p < 0.001; G V vs. G K : t(12) = -8.521, p < 0.001).

. G N C N G N C K G K C N G N C V G V C N G V C K G K C V G N C N .
Finally, subjects were asked to choose the best overall feedback condition and which feedback they preferred for the two provided pieces on information. Eight subjects out of thirteen chose G K C V , four chose G V C K , and one chose G K C N . Nine subjects out of thirteen chose kinesthetic as the best feedback modality for providing guidance, while four chose vibrotactile. Seven subjects out of thirteen chose vibrotactile as the best feedback modality for providing cutting force rendering, while six chose kinesthetic.

V. CONCLUSIONS AND DISCUSSION

We presented an innovative haptic-enabled teleoperation system for the insertion of flexible needles in soft tissue. It combines a state-of-the-art needle tracking and steering framework with a new friction estimation algorithm, with the objective of providing the human operator with informative needle direction guidance and cutting force rendering feedback. To convey both pieces of information in a reliable and easy-to-understand way, we design feedback techniques combining grounded kinesthetic feedback and ungrounded/wearable vibrotactile feedback.

We evaluated the proposed approach in a human subjects experiment enrolling 13 subjects. Participants were asked

G N C N (white), G N C K , G K C N , G N C V , G V C N , G V C K , G K C V (black).
N C N (white), G N C K , G K C N , G N C V , G V C N , G V C K , G K C V (black)
. In (b), from left to right, we have the four information provided conditions C V , C K (red) and G V , G K (blue). The red bars indicate statistically significant differences.

to insert a flexible needle in a gelatine phantom until they reached a stiff target object. We tested the performance of the teleoperation system in seven experimental conditions, i.e., different combinations of kinesthetic/vibrotactile feedback to convey needle guidance and/or cutting force, measuring both objective and subjective metrics. Results show that conditions employing both kinesthetic and vibrotactile feedback always performed better than those employing only one of them. While this result might seem obvious (more information equals better performance), it is often believed that providing multiple pieces of haptic sensations at the same time results confusing -this was not our case. The most relevant performance differences in the objective results appeared in the target accuracy, while no statistically significant difference was found for the other objective metrics, i.e., participants were more precise but not faster. Our results are comparable with those of [START_REF] Rossa | Multiactuator haptic feedback on the wrist for needle steering guidance in brachytherapy[END_REF] in terms of insertion success rate and comparable with [START_REF] Abayazid | Experimental evaluation of co-manipulated ultrasoundguided flexible needle steering[END_REF] and [START_REF] Pacchierotti | Teleoperation of steerable flexible needles by combining kinesthetic and vibratory feedback[END_REF] in terms of insertion accuracy, using visual plus vibratory and kinesthetic plus vibratory feedback, respectively. Another important result is the good match between objective and subjective results. In most cases, conditions that showed good performance in terms of targeting accuracy where also well appreciated by the users. Results also showed that needle direction guidance was significantly better followed/used when provided through kinesthetic feedback with respect to vibrotactile feedback. Conversely, vibrotactile feedback seemed to better render the cutting force, although its benefits with respect to kinesthetic feedback are less evident. This result is again shown both in the objective and subjective metrics. This preference of application for the two types of feedback, i.e., kinesthetic is better/preferred for needle guidance while vibrotactile for cutting force rendering, is not so surprising. Indeed, although vibrotactile feedback has been used for navigation, it is well known that it becomes quickly unpleasant if repeated too often [START_REF] Devigne | Power wheelchair navigation assistance using wearable vibrotactile haptics[END_REF]. As also suggested by one subject, this could be the main reason for a degraded performance. On the other hand, the cutting force rendering algorithm only generated vibrations toward the end of the task, when the target object was punctured. This limited activation span and the alarm/buzz nature of the vibrotactile feedback is probably to thank for the increased accuracy performance in the conditions employing it. Finally, we want to highlight a difference in terms of safety between vibrotactile and kinesthetic guidance. Kinesthetic feedback acts directly on the operator's console, making it harder to move away from the preplanned trajectory; on the other hand, vibrotactile feedback provides rich navigation cues while enabling the operators to easily disregard them if necessary, e.g., if the insertion trajectory must be changed suddenly due to an unexpected event.

  Local systemFig.1. Haptic-enabled teleoperation system. The user controls the robotic manipulator and, in turn, the flexible needle through an Omega.6 haptic interface. Navigation guidance and needle tip cutting force are conveyed through a combination of kinesthetic and vibrotactile feedback, provided by the grounded interface and a wearable vibrotactile device, respectively. Finally, the ultrasound (US) images and a reconstruction of the needle pose are shown on a screen posed in front of the user.

Fig. 4 .

 4 Fig. 4. Friction force estimation. (a) Time vs. force applied on the needle f z s (blue) while applying alternated step velocities v i = ±0.005 m/s (orange) along the needle axis direction. (b) Time vs. force applied on the needle for different velocity magnitudes v i ∈ v during needle extraction (magnification of the red square sector above). (c) Friction per unit length f f for different insertion velocities v.

Fig. 5 .

 5 Fig. 5. Example of needle direction guidance with vibrotactile feedback, G V .Vibrotactile sensations guide the user toward minimizing d tip , i.e., the distance between the needle tip x tip and the projection of the target onto the plane perpendicular to the needle P ⊥ , i.e., x P targ . In this example, as t d < d n,x tip < tu, t d < d n,z tip < tu, -d n,x tip < t d , and -d n,z tip < t d , we calculate m(1), m(2) using the middle term in (5) and m(3), m(4) using the first term.

Fig. 6 .

 6 Fig. 6. Visual feedback showing US images and the needle reconstruction.

Fig. 7 .

 7 Fig. 7. Cutting force f z t and tip insertion velocities vn = [ v x n , v y n , v z n ] in a representative insertion in modality G K C V . The blue-shaded area represents the time when the armband was vibrating.

  Targeting accuracy (cm).

  Insertion success rate (%).

  Insertion speed (mm/s).

Fig. 8 .

 8 Fig. 8. Human-subjects study: objective results. Mean and standard deviation of the mean of (a) targeting accuracy, (b) insertion success rate, (c) insertion time, and (d) insertion speed for the seven conditions -from left to right in the figures -G N C N (white), G N C K , G K C N , G N C V , G V C N , G V C K , G K C V (black).The red bars in (a) indicate statistically significant differences.

Fig. 9 .

 9 Fig. 8. Human-subjects study: objective results. Mean and standard deviation of the mean of (a) targeting accuracy, (b) insertion success rate, (c) insertion time, and (d) insertion speed for the seven conditions -from left to right in the figures -G N C N (white), G N C K , G K C N , G N C V , G V C N , G V C K , G K C V (black).The red bars in (a) indicate statistically significant differences.

* M. Aggravi and D. A. L. Estima contributed equally to this work. This work was supported by funds from The Netherlands Organization for Scientific Research (Innovational Research Incentives Scheme-VIDI: SAMURAI project no. 14855).


Hamid Reza Vanaei
email: hamidreza.vanaei@ensam.eu

Mohammadali Shirinbayan

Sidonie Fernandes Costa

Fernando Moura Duarte

José António Covas

Michael Deligant

Sofiane Khelladi

Abbas Tcharkhtchi

| Fernando

Moura Duarte

Orcid Hamid

Reza Vanaei

Experimental study of PLA thermal behavior during fused filament fabrication

à la diffusion de documents 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.

ment of the filaments during the deposition sequence. This study reports experimental measurements of filament temperature during deposition. These temperature profiles were compared to the predictions made by a previously developed model. The two sets of data showed good agreement, particularly concerning the occurrence of reheating peaks when new filaments are deposited onto previously deposited ones. The developed experimental technique is shown to demonstrate its sensitivity to changing operating conditions, namely platform temperature and deposition velocity. The data generated can be valuable to predict more accurately the bond quality achieved in FFF parts.

K E Y W O R D S

manufacturing, thermal properties, thermoplastics

| INTRODUCTION

Additive manufacturing (AM) denotes a group of innovative technologies that enable the rapid fabrication of three dimensional (3D) physical objects directly from computer-aided design (CAD) data without the use of tooling. Parts with complex geometry that are difficult to produce using traditional manufacturing processes can be obtained by AM. 1,2 A large array of AM techniques is currently used to process thermoplastic polymers, polymer composites, metals, and ceramics. [3][4][START_REF] Saeidlou | Cinétique de Cristallisation, Structure Et Applications Des Stéréocomplexes de PLA[END_REF][START_REF] Santos | [END_REF] Fused filament fabrication (FFF) is extensively used to produce prototypes for applications in, for example, the aerospace, medical, and automotive industries. 7,8 In this process, a thermoplastic polymer is fed into a liquefier that extrudes a filament while moving in successive X-Y planes along the Z direction, to fabricate a 3D part layer-by-layer. [START_REF] Chennakesava | Fused deposition modelinginsights[END_REF] Consequently, as the deposition progresses, the hot filament is deposited onto filaments that were previously deposited and that are now in the process of cooling. This causes their reheating, defining a time during which the interfaces of contacting filaments are above the glass transition temperature (T g ) in the case of amorphous materials, or of the crystallization temperature (T c ) for semicrystalline materials, which is necessary for proper bonding to take place. [START_REF] Hopper | [END_REF]11 Therefore, each filament should be sufficiently hot during deposition, but not too hot, to avert deformation due to gravity and to the weight of the filaments deposited in subsequent layers.

Given the above, the evolution of the temperature profile of the filaments during deposition is a key parameter affecting the bonding quality. 12 Both experimental and theoretical approaches have been proposed to obtain data on the temperature profile of a printed structure. 13 Generally, temperature measurements are limited to a single location. Also, the deposition may be interrupted to fix a thermocouple to the part being fabricated. 14 To circumvent this limitation, infrared (IR) thermography has been used. This approach yields the surface temperature, but it cannot read the temperature at the interface of adjacent layers, both due to the camera measurement accuracy and the small filament dimensions. 15 Ferraris et al. 16 used IR thermography to measure the temperature profile of a vertical wall, but the comparison with theoretical predictions showed poor agreement. Kousiatza et al. 17 applied local measurement of the temperature profile in a specific case on acrylonitrile butadiene styrene (ABS) using K-type thermocouples (d = 250 μm) by pausing the process manually to add them. Although they have concluded that there is a good agreement between experimental and theoretical results, the sudden drop of temperature at the head tip of the extruder clearly shows that there is still a gap between the monitored and numerically derived temperature peak values.

Theoretical efforts to model the temperature profile have focused on 1D or 2D approaches. Sun et al. 18 and Zhang et al. 19 investigated both numerically and experimentally the effect of numerous parameters that exist in the nature of the process. Yardimci and Güçeri 20 developed a 1D numerical model to predict the cooling/bonding in fused deposition of ceramics, considering exclusively convection with the environment. Bellehumeur et al. 21 also developed a 1D model with the same assumption of taking into account the temperature profile along the length of the deposited filament. More recently, Costa et al. 22 developed a computer code that takes into account the heat transfer between adjacent filaments during deposition and predicts temperature and adhesion quality for most 3D-printed parts. The model was shown to be in good agreement with experimental data. 23 Despite the above developments, there is still a lack of practical knowledge on the temperature development of filaments during the deposition stage in FFF. To address this limitation, A. Tcharkhtchi et al. 24 added thermocouples to the build simultaneously with the fabrication without damaging or pausing the process. However, a drop of approximately 50 C was observed on the experimental data. This work presents an improved measurement setup enabling to record the temperature evolution in various locations of the part during deposition, including the interface between two adjacent filaments. The data collected is compared with the predictions provided by the model of Costa et al 22 that have been already developed.

The paper is organized as follows. The experimental procedure is presented in detail and the used heat transfer model is introduced. The temperature evolutions of a single PLA filament and a vertical wall are measured experimentally and the data are compared with the theoretical predictions for validation purposes. Finally, the usefulness of the new experimental method is illustrated by studying the influence of the platform temperature and deposition velocity on the heat transfer during cooling.

| EXPERIMENTS 2.1 | Materials

A commercial PLA filament with a diameter of 1.75 mm (± 0.01 mm) and a density of ρ = 1.24 g/cm 3 have been used. Differential scanning calorimetry (DSC Q10 from TA Company) determined the crystallization and the melting temperature of the material before and after printing. The test was performed by subjecting a 6.8 mg sample from room temperature to 210 C at a heating rate of 5 C/min under nitrogen atmosphere. To measure the main transition temperatures, dynamic mechanical analysis tests were performed (DMA Q800 from TA Company) under tensile mode from 40 to 100 C at a temperature rate of 2 C/min and a frequency of 1 Hz. The rectangular sample with a dimension of 25 × 10 mm 2 was used. For both DSC and DMA characterization, 'TA instrument' software was applied in measurement of the T g , the T c , the melting temperature (T m ) and the enthalpy at these temperatures (ΔH) of both materials.

| 3D printing

The printed parts were manufactured by a desktop 3D printer fitted with a single nozzle (d = 0.4 mm) printing head and a temperature-controlled atmosphere (build platform). The solid model file corresponding to the part illustrated in Figure 1 was designed using the FreeCad software and then exported as stereolithography (STL) format to be loaded into the FlashPrint software, which generates the printing path.

The test case was built considering the values for the processing variables (Table 1) that are commonly used in the desktop 3D printer to ensure a good quality part in terms of bonding between filaments and mechanical strength. 21,[25][26][27][28]

| In situ temperature measurements

In situ temperature measurement methods should be sufficiently precise and quick to track filament cooling and the reheating peaks arising from contact between freshly and previously deposited filaments. In addition, it should be possible to apply the sensor locally without the need to pause the process. A very small K-type thermocouple (d = 80 μm), capable of measuring temperatures between -75 C and 250 C, was used. By taking advantage of this size, it is possible to squeeze the device between two adjacent filaments. The thermocouple is connected to a Datapaq® Tracker Telemetry system (previously employed for the in situ temperature measurement in rotational molding process 29 ) for temperature recording and connection to a polycarbonate (PC). As shown in Figure 2 (setup of the work), when the print was started by fabrication of the support (by opening the door) thermocouple placed at the location of 5 mm from the starting point of the deposition and then the door closed in order to have a stabilized temperature of the environment.

| MODELLING OF FILAMENT COOLING DURING DEPOSITION

Heat transfer during deposition is complex, with contributions from radiation, convection, and conduction. However, it has been demonstrated that (a) the heat losses by convection with the environment (b) the thermal contacts with the support and with adjacent filaments are the main contributors to the filament temperature evolution. 22 A computer code that has been already developed, was applied assuming the gradual deposition of small axial filament segments, an analytical solution for the energy equation 30 whilst updating the local thermal conditions, and a healing criterion proposed by Yang and Pitchumani. 31 This gave rise to a useful tool that allows us to predict the temperature evolution and the degree of bonding between filaments for 3D parts including the usage of two distinct materials (e.g., the material of the part plus support material).

Figure 3 shows the evolution of the temperature of the filaments at specific instants upon building the first ten layers of a vertical wall. As in the experiments

| Validation of the measurement methodology

For validation purposes, the temperature evolution of PLA during the deposition of a single filament and of a vertical wall as measured using the applied methodology was compared with the predictions of the thermal model.

| Deposition of a single filament

The temperature evolution during the deposition of a single filament is presented in Figure 5 for the first deposited filament. Both experimental measurements and theoretical curves for two values of the heat transfer coefficient are shown. Practice revealed that opening/closing the door of the environmental chamber to add/remove thermocouples would slightly disturb the platform temperature. This is why an interval in the range 50-60 C (approximated using thermocouples to see the temperature variation of the platform) and not the value of 60 C (see Table 1) is shown in the mentioned figure (Figure 5). The value of 70 W/m 2 C for the heat transfer coefficient (h conv ) is commonly used (e.g., 32 ). A value of 88 W/m 2 C is obtained when using the Churchill correlation for the cooling down of a cylinder by natural convection 33 :

h conv = Nu d :k d : ð1Þ
where d is the diameter (m), k is the thermal conductivity (W/m C) and Nu d is the Nusselt number defined by:

Nu d = 0:60 + 0:387 Ra 1=6 d 1 + 0:559 Pr À Á 9=16 h i 8=27 8 > < > : 9 > = > ; 2 : ð2Þ
where the Rayleigh Ra d number and P r are expressed as:

Ra d = Gr d Pr Pr = v k α ( : ð3Þ
In the above expressions, v k is the kinematic viscosity (m 2 /s), α is the thermal diffusivity (m 2 /s), and Gr d is the Grashof number:

Gr d = gβ T s -T E ð Þ d 3 v 2 k : ð4Þ
Here g is the gravity acceleration (9.8 m/s 2 ), β is the volumetric thermal expansion coefficient, T s is the cylinder temperature ( C) and T E is the environment temperature ( C).

Figure 5 shows a good agreement between the theoretical and the experimental values. The difference between the two sets of data occurs mostly between 3 and 8 s, when the predicted cooling rate is higher than the one measured. This is probably due to the fact that the theoretical model does not consider the change in state and crystalline growth, and thus forecasts faster cooling. reheating of the filament due to the deposition of filaments 2 to 4). The crystallization and platform temperatures are also identified. For computational purposes, it is important to define the thermal contact conductance (h) between adjacent filaments. This is difficult, as it depends on pressure, surface roughness, and other conditions that are difficult to quantify. [START_REF] Holman | Transferência de Calor. S ão Paul o[END_REF] Apparently, there are no theoretical or empirical correlations providing an exact value for h. Using a value of h = 800 W/m 2 C, the magnitude of the experimental and theoretical reheating peaks became virtually coincident.

Regardless of this approximation, the onset, relative magnitude, and breadth of the various temperature peaks are similarly captured by the two approaches. As expected, the peaks become gradually smaller with time, as the new filament being deposited is separated from the first filament by more filaments. As before, the predictions seem to overestimate the cooling rate, as no phase change and crystallization were built in the model. On the other hand, a delay in receiving the experimental data can exist and contribute to the differences.

| Case studies

This section demonstrates the usefulness of the proposed in situ temperature measurement technique, by studying the effect of the platform temperature and deposition velocity on the heat transfer during cooling.

Figure 7 shows the temperature evolution of three different temperature of the platform, at location distance 5 mm from the edge. As expected, the lower the platform temperature the faster the cooling. As for the reheating peaks, they have identical onsets, but the magnitude tends to decreases with increasing platform temperature. When the platform temperature is set to 100 C, the filament being monitored reheats repeatedly above its T c , favoring bonding.

When the deposition velocity increases (Figure 8), the rate of cooling decreases. Also, and as expected, the onset of the peaks occurs at different times and their breadth is also altered. Finally, the magnitude of the peaks for the lowest deposition velocity is higher, probably due to the higher difference between the temperatures of adjacent filaments.

| CONCLUSIONS

This work presents a 'localized measurement' setup enabling to record temperature profile of two adjacent filaments (or a sequence of deposition) during the FFF process. The main characteristic of the proposed method is the accurate measurements of polymer material temperature using 80 μm thermocouples. It has been shown that the experimental results are in good agreement with the utilized numerical method. It correctly measures the temperature profile of a single filament and filaments embedded in multilayer parts by taking into account the effect of the variant temperature of the environment and thermal contact conductance between adjacent filaments. A parametric study on the effect of the bed temperature and the printing speed indicated the influence on reheating of the previously deposited filament during additional layer deposition.

These results are important for the understanding of filament bonding. The main outcome of this study is to apply the incoming results for prediction of bonding and consequently to optimize the strength of successively deposited layers. Future work will focus on the improvement of the proposed setup to show its capacity in complex geometry. Accordingly, it could help the optimization of bonding quality by correlating the thermal and mechanical characteristic results.

F 4 |

 4 I G U R E 1 Schematic of the test case representing the deposition mechanism of each layer T A B L E 1 Process parameters used for printing U R E 2 (a) Setup and b) schematic used for the in situ measurement of temperature during the deposition stage in FFF [Color figure can be viewed at wileyonlinelibrary.com] reported in this work, the liquefier deposits one filament, then stops and returns to the initial point to deposit the next filament. Under these printing conditions, when a new filament is deposited, the previous one has already significantly cooled down. Nevertheless, the deposition of a new hot filament prompts the reheating of filaments of the previous layers (as seen at 32.5 and 41.5 s), thus demonstrating the importance of considering the thermal contacts in the calculations. RESULTS AND DISCUSSION 4.1 | Polymer characterization

Figure 4 (

 4 Figure 4(a) shows the DSC traces for the PLA filament and the printed part. The T g , the T c , the T m , and the

F I G U R E 3

 3 Temperatures of the ten first layers of the vertical wall at some instants of the deposition process [Color figure can be viewed at wileyonlinelibrary.com]

F

  Abbreviations: DMA, dynamic mechanical analysis; DSC, differential scanning calorimetry.

Figure 6

 6 Figure6depicts the evolution of the temperature of the first filament (at a location distant 5 mm from the edge) during the building of a vertical wall consisting of single filaments deposited on top of each other. It is clear that cooling of this first filament is significantly affected by the successive deposition of younger filaments, which may cause important reheating. The numbers identify regions of the data (1: cooling of the first filament; 2-4: reheating of the filament due to the deposition of filaments 2 to 4). The crystallization and platform temperatures are also identified.For computational purposes, it is important to define the thermal contact conductance (h) between adjacent filaments. This is difficult, as it depends on pressure, surface roughness, and other conditions that are difficult to quantify.[START_REF] Holman | Transferência de Calor. S ão Paul o[END_REF] Apparently, there are no theoretical or empirical correlations providing an exact value for h. Using a value of h = 800 W/m 2 C, the magnitude of the experimental and theoretical reheating peaks became virtually coincident.Regardless of this approximation, the onset, relative magnitude, and breadth of the various temperature peaks are similarly captured by the two approaches. As expected, the peaks become gradually smaller with time, as the new filament being deposited is separated from the first filament by more filaments. As before, the predictions seem to overestimate the cooling rate, as no phase change and crystallization were built in the model. On the other hand, a delay in receiving the experimental data can exist and contribute to the differences.

F

  I G U R E 6 Experimental (± 2 C) and theoretical temperature evolution of the first filament (at x = 5 mm) during the deposition of a vertical wall consisting of single filaments deposited on top of each other. The numbers identify regions of the data: 1-cooling of the first filament, 2-4 reheating of the filament due to the deposition of filaments 2-4 [Color figure can be viewed at wileyonlinelibrary.com] F I G U R E 7 Temperature profile (± 2 C) of vertical wall at x = 5 mm from the start of deposition at different platform temperatures [Color figure can be viewed at wileyonlinelibrary.com] F I G U R E 8 Temperature profile (± 2 C) of vertical wall at (x = 5 mm) from the start of deposition at different print speed [Color figure can be viewed at wileyonlinelibrary.com]