Samuel Veilleux
email: samuel.veilleux.4@ulaval.ca

Najat Bhiry

Armelle Decaulne

Talus slope characterization in Tasiapik Valley (subarctic Québec): Evidence of past and present slope processes

Keywords:

Topographic, granulometric, morphometric, petrographic and vegetation surveys were conducted on the slopes of Tasiapik Valley, near Umiujaq (Nunavik), to document mass wasting processes and their geomorphological impact. Talus slopes, widespread at the foot of the steep rockwalls of Tasiapik Valley, are an important landscape feature in the area. The lithology of the slope deposits attest their local origin, namely the result of rockfalls coming from the adjacent wall. Locally, poor vegetation covering the clasts exhibits recently fallen debris; elsewhere, dense shrub cover has colonized the slopes demonstrating the low activity nowadays. On-going periglacial processes have led to extensive dismantling of the rockface, enabling for debris supply. Following the last deglaciation, paraglacial processes have potentially favoured slope instabilities. The use of automatic cameras during the winter 2017-2018 resulted in the observation of many snow-avalanche events; however few rockfall events have been observed. Spring snow avalanches have carried rock debris to the talus at the foot of the slope; snow also enabled debris redistribution on the slopes.

Introduction

Northern landscapes have undergone many changes since their deglaciation. In particular, paraglacial conditions (e.g. glacioisostatic rebound and rockface dismantlement) induced talus slope formation by supplying debris through the pressure release on rock fractures and freeze-thaw processes [START_REF] Ballantyne | Paraglacial Slope Adjustment and Resedimenfation following Recent Glacier Retreat, Fåbergstølsdalen, Norway[END_REF][START_REF] Matsuoka | Rockfall activity from an alpine cliff during thawing periods[END_REF][START_REF] Ballantyne | Paraglacial geomorphology[END_REF][START_REF] Matsuoka | Frost weathering and rockwall erosion in the southeastern Swiss Alps: long-term (1994-2006) observations[END_REF].

Nunavik is part of the low-Arctic region of eastern Canada and its landscape consists of low hills, basins and plateaus. The few studies that have been conducted in this vast region have demonstrated the occurrence of slope processes on slopes less than 100 m high [START_REF] Belzile | Les versants rocheux périglaciaires à la presqu'île des Manitounouc Kuujjarapik[END_REF][START_REF] Bégin | Analyse dendrochronologique d'un glissement de terrain de la région du Lac à l'Eau Claire (Québec nordique)[END_REF][START_REF] St-Cyr | Formation et évolution des versants rocheux des îles centrales du lac à l'Eau Claire[END_REF][START_REF] Marion | The Holocene development of a debris slope in subarctic Québec, Canada[END_REF][START_REF] Germain | The vulnerability of northern cities to weather-related hazards: case studies from the Province of Quebec, eastern Canada[END_REF][START_REF] Germain | Snow avalanche hazard assessment and risk management in northern Quebec, eastern Canada[END_REF]. Recent studies [START_REF] Decaulne | Geomorphic evidence of Holocene slope dynamics on the Canadian shield-a study from Lac à l'Eau-Claire, western Nunavik[END_REF][START_REF] Bhiry | Development of a subarctic peatland linked to slope dynamics at Lac Wiyâshâkimî (Nunavik, Canada)[END_REF] conducted at Wiyâshâkimî Lake in Tursujuq National Park (Nunavik) showed that talus slope formation started after deglaciation at about 4600 BP, and that slope processes are still active today. Some of the villages in Nunavik (Salluit, Kangiqsujuaq and Kangiqsualujjuaq) are located within glacial valleys with prominent slopes, while other villages (Umiujaq) are situated near high cuesta relief (∼230 m). Accordingly, it is crucial to document slope dynamics and to evaluate associated risks on the local population, visitors and infrastructures. For instance, in Kangiqsualujjuaq (northeastern Nunavik), a dreadful snow avalanche hit the gymnasium of Satuumavik school during the 1999 New Year's Eve celebrations, causing the death of 9 people and injuring 25 [START_REF] Bérubé | Rapport d'enquête publique sur les circonstances des décès survenus à Kangiqsualujjuaq, Nouveau-Québec le premier Janvier[END_REF][START_REF] Lied | Avalanche Hazard Assessment in Nunavik and on Cotenord[END_REF][START_REF] Germain | Snow avalanche hazard assessment and risk management in northern Quebec, eastern Canada[END_REF]. However, no extensive research has been conducted on slope processes in the Umiujaq area (including snow avalanches, landslide and rockfalls), their triggering factors, their occurrence and their runout distance. General conditions are conducive for bedrock dismantling and mass wasting, even with limited slope heights, but additional knowledge about slope processes is still required.

The main objective of the study is to document landforms organisation built by gravitational processes in Tasiapik Valley and their contribution to talus development based on geomorphological surveys. This study discusses slope evolution during the Holocene, from the retreat of the Laurentide Ice Sheet in the area to the present-day, highlighting the potential risk at the valley bottom. 1.1. Regional setting Tasiapik Valley (56 • 33 N, 76 • 28 W) is located 5 km east of the Inuit village of Umiujaq, on the east coast of Hudson Bay in Nunavik, Québec (Fig. 1a). It is approximately 4.5 km long and 1.5 km wide, following a northwest-southeast orientation. At the southeastern end of the valley lies Tasiujaq Lake (formerly named Guillaume-Delisle Lake or Richmond Gulf), a 691 km 2 brackish water body connected with the Hudson Bay by a narrow cataclinal channel called Le Goulet (ARK, 2007). The lake is part of Tursujuq National Park, created in 2013.

The regional geology is characterized by a Paleoproterozoic volcano-sedimentary sequence lying unconformably on the Precambrian shield (Fig. 1b). The volcano-sedimentary sequence includes limestone, quartz arenite, dolomite and sandstone strata (Qingaaluk Formation) underlying a thick (∼15 m) basalt layer (Nastapoka Group) dipping westward [START_REF] Stockwell | Géologie du Bouclier canadien[END_REF][START_REF] Chandler | Tectonics of the Richmond Gulf Area, Northern Quebec -a hypothesis[END_REF][START_REF] Chandler | The Early Proterozoic Richmond Gulf Graben, East Coast of Hudson Bay[END_REF]Eaton and Derbyshire, 2010). This asymmetrical monoclinal relief (cuesta) consists of a gentle western slope and a steep eastern slope and extends over 650 km along the east coast of Hudson Bay (Dionne, 1976;Guimont and Laverdière, 1980). Tasiapik Valley lies at the frontslope of the cuesta on its southwestern side, whereas the northeastern side consists of a residual butte called Umiujaaluk Hill.

Quaternary deposits on the east coast of Hudson Bay are the result of a succession of sedimentary environments following the retreat of the Laurentide Ice Sheet at about 8200 cal. BP [START_REF] Hillaire-Marcel | La déglaciation et le relèvement isostatique sur la côte est de la baie d'Hudson[END_REF][START_REF] Allard | La déglaciation d'une partie du versant hudsonien québécois: bassins des rivières Nastapoca, Sheldrake et à l'Eau Claire[END_REF][START_REF] Lavoie | Deglaciation landforms and C-14 chronology of the Lac Guillaume-Delisle area, eastern Hudson Bay: a report on field evidence[END_REF] (Fig. 1b). Lowlands below 271 m a.s.l, the altitudinal limit of the postglacial Tyrrell Sea in Tasiujaq Lake area [START_REF] Fraser | Morphology and facies architecture of a falling sea level strandplain, Umiujaq, Hudson Bay, Canada[END_REF][START_REF] Lavoie | Deglaciation landforms and C-14 chronology of the Lac Guillaume-Delisle area, eastern Hudson Bay: a report on field evidence[END_REF], are covered by deep-water and shallow-water marine sediments, and littoral deposits (raised beaches) that were formed during stages of rapid glacio-isostatic uplift (Hillaire-Marcel, 1976).

A glaciomarine fan complex lies in the upstream part of the valley. It consists of fluvioglacial material that was deposited during a stillstand of the ice margin around 8000 cal. BP (Lajeunesse and Allard, 2003b).

The study area has a cold subarctic climate and it is located in the discontinuous permafrost zone [START_REF] Allard | Nunavik and Nunatsiavut: From Science to Policy: An Integrated Regional Impact Study (IRIS) of Climate Change and Modernization[END_REF]. Mean annual air temperature recorded between 2013 and 2017 varies between -5.6 and -4.2 • C, with maxima of 23 • C and minima of -36 • C [START_REF] Fortier | Groundwater monitoring network from the Umiujaq region in Nunavik[END_REF]. Mean annual precipitation is approximately 500 mm, with 40% falling as snow [START_REF] Ménard | Monitoring of ground surface temperatures in various biophysical micro-environments near Umiujaq, eastern Hudson Bay, Canada[END_REF]. The Umiujaq area is located at the edge of the shrub and forest tundra zones; low shrubs, ericaceous plants and lichens cover the upstream part of Tasiapik Valley, while dense forest cover occupies the downstream part [START_REF] Payette | The forest tundra and present tree-lines of the northern Québec-Labrador peninsula[END_REF]. Shrub cover has expanded significantly (shrubification) during the 20 th century [START_REF] Ménard | Monitoring of ground surface temperatures in various biophysical micro-environments near Umiujaq, eastern Hudson Bay, Canada[END_REF][START_REF] Provencher-Nolet | Quantification des changements récents à l'écotone forêt-toundra à partir de l'analyse numérique de photographies aériennes[END_REF][START_REF] Pelletier | Ecosystem changes across a gradient of permafrost degradation in subarctic Québec (Tasiapik Valley, Nunavik, Canada)[END_REF].

The SW side of Tasiapik Valley has a near-vertical rockwall. It is approximately 50 m high in the upstream part of the valley, but it increases to 230 m near Tasiujaq Lake in the downstream part. Slope deposits lie at the base of the escarpment, connecting the rockwall to the valley floor, but have also accumulated on basaltic rocky outcrops in the uppermost part of the rockwall (Fig. 1c). The NE side features a step-like topography, with slope deposits either located at the base of the slope or perched on basaltic and sedimentary rocky outcrops. A gravel road connecting Umiujaq to Tasiujaq Lake follows the cuesta frontslope on the SW side.

Methods

For this study, 18 talus slopes were investigated, on the SW and NE sides of Tasiapik Valley (Fig. 1c). Data were collected over four field campaigns during the summers of 2016 (August), 2017 (August) and 2018 (June and August). Several slope deposits were identified by satellite imagery prior to initiating fieldwork. Two sets of orthophotos from Québec's Ministère des Ressources naturelles et de la Faune (MRNF) were used, one dating from 2004 (scale 1/10,000, 25 cm resolution) and the other from 2010 (scale 1/10,000, 15 cm resolution).

Topographic surveys were conducted along 18 longitudinal transects using a Leica DGPS (Differential Global Positioning System). Waypoints were recorded from the apex of the slope deposits to their base, perpendicular to the rock face. Data were processed in ArcGIS and Excel to produce topographic profiles, revealing the microtopographic features such as inflection and texture in accurate details. The estimation of the stage of evolution is carried out from the Ho/Hi ratio, where Ho corresponds to the height of the talus slope and Hi to the total height of the slope including the rockwall [START_REF] Francou | L'éboulisation en haute-montagne -Andes et Alpes -, six contributions à l'étude du système corniche-éboulis en système périglaciaire[END_REF][START_REF] Sellier | Évolution comparée de versants quartzitiques des Highlands d'Écosse et de Norvège centrale (Rates of quartzitic slopes evolution in the scottish Highlands and central Norway)[END_REF]. The ratio gives an overview of the exhaustion of the remaining rockwall (debris source) in concomitance with the talus slope formation on a longer timescale, namely since the last deglaciation. For example, a ratio approaching 1 indicates an advanced slope development stage due to the low height of the residual wall compared to the height of the talus slope.

Topographic data and satellite imagery were used to measure rockfall runout distances (horizontal travel distance) calculated from the source area to the farthest slope debris. In addition, the reach angle, calculated from the source-area to the farthest slope debris, and the shadow angle, calculated from the apex of the talus to the farthest slope debris, were documented to provide information about the extent of slope processes in the area.

On 12 of the 18 longitudinal transects, granulometric and petrographic surveys were conducted by sampling 25 randomly selected rock fragments at intervals of 10-15 m along the transects. Debris were measured along three axes: length (a-axis), width (b-axis) and thickness (c-axis). Measurements were compiled in Excel and then analyzed to produce descriptive statistics. Morphometric indices were also calculated from these measures (Pérez, 1989, Hétu and[START_REF] Hétu | Effects of environmental change on scree slope development throughout the postglacial period in the Chic-Choc Mountains in the northern Gaspé Peninsula, Québec[END_REF]. The flattening index (Fi) is calculated as follows:

F i = a + b 2c (1)
where a corresponds to the length, b to the width and c to the thickness of the fragment [START_REF] Cailleux | L'indice d'émoussé des grains de sable et grès[END_REF]. A high Fi value indicates that the debris has a flatter shape. The elongation index (Li) is calculated as follows:

L i = a b (2)
where a and b correspond to the length and width of the fragment [START_REF] Schneiderhöhn | Eine vergleichende Studie über Methoden zur quantitativen Bestimmung von Abrundung und Form an Sandkörnern (Im Hinblick auf die Verwendbarkeit an Dünnschliffen[END_REF]. A high Li value indicates that the debris tends to be elongated. Finally, the sphericity index (Si) is calculated as follows:

S i = bc a 2 1 3 (3)
where a, b and c correspond to the length, width and thickness of the fragment [START_REF] Krumbein | Measurement and geological significance of shape and roundness of sedimentary particles[END_REF]. A value Si approaching 1 indicates that the debris has a more massive shape, spherical in the case of a rounded fragment and cubic for an angular fragment. These indices document the falling behavior of clasts, since spherical debris are prone to rolling, while elongated flat debris are more likely to slide. Petrographic surveys provide lithological data for the measured fragments. Their origin, either local (associated with the local slope development) or exogenous (generally from glacial transport and deposition), is closely related to their lithology, thus their general shape, and their position on the slope. The edges of the debris were characterized, with a view to determining their origin: an angular fragment has undergone very little erosion, indicating a short transportation distance/local source, while debris transported by glaciers or reworked in the Tyrrell Sea has a pronounced rounded shape.

Vegetation cover was described at each sampling station in order to assess recent and current process activity. Hierarchical values were attributed to each station based on the type of vegetation and the estimated percentage of coverage on the debris, providing relative age-estimates: 1) Fresh debris: no lichen species observed on the clast; 2) Recent debris: some lichen species observed on the clast; 3) Medium-aged debris: several lichen species partially cover the clast; 4) Old-aged debris: several species of lichens and mosses partially cover the clast; 5) Very old-aged debris: several species of lichens and mosses totally cover the clast; potentially also covered with low shrubs.

Vegetation classification values and Ho/Hi ratio values were used to estimate the stage of slope development. The addition of these two values gives an overview of the slope evolution from both short term (vegetation) and long term (Ho/Hi) perspectives. Values of 1 to 5 were assigned to each longitudinal profile according to their Ho/Hi ratio, following the Jenks natural breaks classification method [START_REF] Jenks | The data model concept in statistical mapping[END_REF]; a value approaching 1 indicates a low Ho/Hi ratio, thus a younger development stage. The vegetation values (i.e. the lowest -freshest -value per profile, ranked 1 to 5 according to the relative age estimate described above) were added to provide an overall development score. In addition, the age of shrubs at the bottom of talus along the SW-07, SW-08 and SW-09 profiles was determined using dendrochronology on 11 black spruce samples (Picea mariana).

Finally, in order to monitor slope movements on a shorter time scale, three Reconyx PC800 automatic time-lapse cameras were installed on the SW side of the valley in August 2017. One of the cameras (TAS 1) is located on the cuesta edge, above the rockwall and talus along the SW-07 and SW-08 profiles (Fig. 2). The latter is covered by a second camera (TAS 2) that is located 300 m away from the rockwall. A third camera (TAS 3) is located further upstream near the base of the talus along the SW-06 profile. Approximately 14,000 photos were taken over a one-year period from August 2017 to August 2018 in the valley. Photos were taken during daytime at one-hour intervals until June 2018, then at 15 or 30 min intervals (depending on the location) until August 2018.

Results

Topography of slope deposits

The SW profiles show a steeper slope gradient, with a mean angle of 25.3 • and a median angle of 26.3 • , while the NE profiles have a mean angle of 21.2 • and a median angle of 17.5 • (Table 1; Fig. 3). A slight or distinct concavity of the talus slope is apparent on six of the nine SW profiles, while the SW-01 and SW-06 profiles are virtually linear. The SW-09 profile shows a more complex shape (linear proximal part and chaotic distal part). On the NE side, there is no distinct concave profile, yet seven of the nine profiles show either a slightly concave or a linear shape, while the NE-04 and NE-07 profiles respectively show a complex and a convex shape. Four profiles (SW-07, SW-08, SW-09 and NE-09) exhibit a strong basal concavity (Fig. 3).

Relative dating of slopes

The mean Ho/Hi index on the SW (0.36) and NE sides (0.25) indicates that the remaining rockwall is generally higher than talus slopes (index < 0.5) (Table 2). However, the step-like topography on the NE side, compared to the near-vertical rockwall on the SW side, could mean that the NE side has reached a more advanced stage of development. Talus slopes near the southern margin of both sides (along the SW-07, SW-08, SW-09 and NE-09 profiles) are located under high vertical rockwalls with Ho/Hi index below 0.2, thus indicating a younger development stage. The SW-04 and SW-05 profiles indicate an older stage than the other SW profiles, with respective index values of 0.56 and 0.62. Examination of the vegetation covering the surficial debris on the slope deposits revealed the presence and position of few fresh deposits. Most of the debris had a clear vegetation cover on the SW side (Fig. 4). In the apical parts of the talus slopes, various lichens and/or mosses are abundant, while a discontinuous thin strip of herbaceous plants and low shrubs is located at the edge of the rockwall undisturbed by present-day slope activity. Distal parts also feature abundant lichens and mosses on most of the clasts along with thick mosses covering the slope deposits and low shrubs; this trend is especially evident on the SW-07 and SW-09 profiles. A similar trend was observed on the NE profiles, as most sampling stations in the apical part of the talus slopes show medium to old-age status, while sampling stations in the distal parts indicate older-age status. However, debris along the SW-08, NE-06 and NE-09 profiles appear to be more recent, with little overall coverage and the presence of few lichens. Some fresh debris were scattered along most of the profiles (Table 3).

By adding up the Ho/Hi index values and the vegetation classification values, we can estimate the developmental stage of the slopes. As shown on Fig. 5, the SW-08 and NE-09 profiles seem to be at the youngest developmental stage among all the talus slopes. The Ho/Hi index and the vegetation classification values are consistent for some profiles, showing a concomitance for both parameters. For example, the talus slope along NE-09 profile has a low Ho/Hi index (0.183) and there is very poor lichen cover on the debris. However, the two parameters proved to be contradictory for some talus slopes, particularly for the SW-09 profiles, due to the overly high rockwall (increased debris supply potential) and the presence of well-developed vegetation (limited debris supply on the talus slope). Both of these findings indicate that the debris supply is sporadic.

Source and morphometry of slope deposits

Three classes of debris were identified along the investigated profiles: 1) basalt, 2) sedimentary rocks (comprising dolomite, limestone, quartz arenite and sandstone), and 3) granitic gneiss (Fig. 6).

On the SW side, sedimentary rock debris comprised 64.5% of the sampled clasts, whereas basalt debris accounted for 35.4% and gneiss for 0.07%. Sedimentary rock debris represented a larger proportion on the NE side, accounting for 89.4% of all clasts, while basalt and gneiss accounted for 6.4% and 4.2%. Assuming the top basalt layer is ∼15 m thick throughout the valley; those values coincide with the large proportion of sedimentary rock strata available for debris supply on the exposed rockwall. Sedimentary rock strata account for 90% (∼140 m) of the rockwall (∼155 m) above the SW-07, SW-08 and SW-09 profiles. In the upstream part of the SW side, sedimentary rock strata account for 57% (∼20 m) of the rockwall (∼35 m). On the NE side, the top basalt layer has considerably receded on most of the investigated slopes, revealing rocky outcrops composed of sedimentary rock strata. However, the rockwall above the NE-06 profile is mainly composed of basalt (68% of total height).

By comparing the respective proportions of each lithology at sampling stations along the profiles, many of the slope deposits show consistent ratios of sedimentary rocks and/or basalt debris from the apex to the base of the slope. On the SW side, the proportion of basalt debris varies between 80% and 84% throughout the sampling stations on the SW-02 profile and between 14% and 26% on the SW-07 profile. On the NE side, the proportion of sedimentary rock debris varies between 80% and 88% on the NE-02 profile, while the NE-08 and NE-09 profiles show no difference as sedimentary rock debris compose 100% of the talus. However, some of the slope deposits show an increasing proportion of basalt material toward the foot of the talus. For example, along the SW-08 profile, the basalt debris percentage increases from 17% at the apex to 37% at the bottom of the slope; on the SW-09 profile, it increases significantly from 12% to 100%. Finally, on the NE-06 profile, the basalt debris are only located at the bottom of the talus, whereas the sedimentary rock debris comprises the entire apical part. Clasts were measured to document the size and morphometry of the slope deposits, and their distribution along the deposits (Fig. 7). On the SW side, the debris (a-axis) range in size from 30 to 154 cm on average (Table 3). However, standard deviation values show important disparities between the measured clasts, as the largest clasts range from 200 to 1900 cm. In addition, mean flatness index values range from 2.36 to 4.22. Lengthening and sphericity indices do not vary much, ranging from 1.48 to 1.86 and from 0.54 to 0.65. On the NE side, sizes range from 83 to 178 cm on average, while standard deviation values show important disparities between the along all profiles, the clasts tend to have a flatter and more elongated shape in the downstream part of the valley on the SW side (SW-07, SW-08 and SW-09), and along the NE-06 and NE-09 profiles. This morphology is more associated with sedimentary rock debris due to the dismantling of the thin sedimentary layers which represents the majority (84.5%) of the sampled debris along these profiles (given the larger proportion of sedimentary rock within the rockwall). For most of the profiles, the mean and median a-axis values tend to increase toward the base of the slopes. As for the morphometry indices, 50% of the profiles show an increasing trend for flatness and sphericity values toward the base of the slope, while the lengthening index of clasts increases from the apex to the base for 75% of the profiles. Debris from every sampling station on the SW side were analyzed and described as being either angular or subangular, which indicates their local slope-related provenance, whereas the NE side shows a greater diversity of debris shape. For the NE-02 and NE-04 profiles, rounded/sub-rounded debris are abundant at the base and angular debris are found near the apex of the talus slopes, meaning that slope-related debris have accumulated in the proximal part and mixed with rounded heavily reworked debris toward the distal part.

Debris runout

Several scattered clasts are located in the distal parts of the slopes. These are mostly located on the SW side and some of them were deposited only a few meters from the road. Their angular shape and lithology (mostly basalt) indicate their local sloperelated provenance, by contrast with rounded fluvioglacial/glacial debris.

According to [START_REF] Corominas | Integrated landslide susceptibility analysis and hazard assessment in the principality of Andorra[END_REF], the maximum reach angle for small-scale (1-10 m 3 ) rockfalls on an unobstructed path is 48 • . Therefore, material falling from the uppermost source area in the rockwall (basalt layer) would theoretically be transported no further than the talus at the bottom of the slope, where the terrain levels out and the shrub vegetation is often denser (Fig. 8). This perfectly matches the very large boulder accumulation on the SW-09 profile. However, the reach angle for the farthest basalt debris, located beyond the talus slopes, ranges from 24 • to 40 • . These lower reach angles show a greater horizontal displacement of fallen debris, resulting from either a large-scale rockfall (10-100 m 3 for a 40 • reach angle; 100-1000 m 3 for a 33 • reach angle; >1000 m 3 for a 26 • reach angle) [START_REF] Corominas | Integrated landslide susceptibility analysis and hazard assessment in the principality of Andorra[END_REF] or an external process. For example, there are a dozen large basalt boulders (a-axis > 100 cm) located in the distal part of the SW-07 profile at a 30 • angle. However, the volume is not sufficient for such a displacement to result from a rockfall. The shadow angle ranges from 10 • to 25 • , with the lowest angles measured on the SW-07 (10 • ) and SW-09 (17 • ) profiles. These angles are lower than the minimum rockfall shadow angle ranging from 22 • to 30 • [START_REF] Rapp | Recent development of mountain slopes in Kärkevagge and surroundings, northern Scandinavia[END_REF][START_REF] Govi | Photo-interpretation and mapping of the landslides triggered by the Friuli earthquake (1976)[END_REF][START_REF] Lied | Rockfall problems in Norway[END_REF][START_REF] Hungr | Engineering evaluation of fragmental rockfall hazard[END_REF][START_REF] Evans | The assessment of rockfall hazard at the base of talus slopes[END_REF]. Therefore, it can be assumed that slope debris were deposited beyond the rockfall runout zone by another process. The hypothetical travel distance for the farthest debris ranges from 105 to 318 m, with the longest distances measured in the downstream part of the SW side.

In June 2018, numerous dirty snow-avalanche deposits were observed on the SW side, but their terminus rarely exceeded the talus slope, thus supporting the assumption that snow-avalanche runout is generally limited to the foot of the talus. However, longer runout was observed for clean snow avalanches occurring in February through April. Extreme snow avalanche runout cannot be documented at this stage in the valley.

Tree-ring counting on 11 black spruce (Picea mariana) samples, dating back to 1900-1930, and located at the foot of the SW-07 and SW-08 profiles, showed clear signs of eccentric growth, such as the formation of reaction wood during their lifespan. As these reaction wood periods are not significantly concomitant from one tree to another, wind or snow cover could hardly be the main controlling factors for these deformations. Instead, slope processes, such as snow avalanches, are likely causes. However, the limited number of samples does not support a precise chronology of events.

Short-term slope dynamics

Analysis of the photographs from the three automatic cameras has documented the slope dynamics from summer 2017 to summer 2018, with the record of one full winter. Snow avalanches were the main gravitational processes observed, occurring from November 2017 to June 2018; snow-avalanche events occurred more frequently from the end of May 2018. During this period, most of the snow-avalanche deposits exhibited a dirty appearance, because rock debris were incorporated in the dense and humid snow. This characteristic implies a flowing motion with frequent contact with the regolith (Fig. 9a). These observations were validated in the field in June 2018, when about 20 wet snow avalanches were observed on the SW side of the valley, from June 11 to June 15.

Along the SW-07 and SW-08 profiles (TAS2), some larger rock debris (∼1 m a-axis) were transported and deposited downslope by snow avalanches. Other debris were observed falling onto the snow-covered talus until early July 2018, concomitant with snowavalanche events. Two significant rockfalls occurred one-hour apart on June 30 2018. Several debris could be observed sliding on the snow covered talus after their fall; they travelled ∼30 m before settling at mid-slope (Fig. 9b). During the same period, only a few small snow avalanches occurred in the area covered by camera TAS3, along the SW-06 profile, and no rockfall were observed. Out of the seasonal presence of snow, no movement was observed on the slopes from the analysis of the photographs.

Based on the observations during 2017-2018, slope movements occur more frequently in the spring. Prior to June 2018, no rocky deposits on the snow cover had been observed following snowavalanches, which occurred sporadically from December 2017 to April 2018.

Discussion

Talus slope formation

Results from the morphometric and petrographic surveys suggest an accumulation of rock debris on the SW and NE sides resulting from successive discrete rockfalls that formed into scree slopes. The vast majority of the sampled debris is from a local source, namely the lithologies exposed on the rockwall, and they exhibit subangular to angular shapes.

Vegetation cover on the talus suggests that the most recent debris supply occurred along the SW-08, NE-06 and NE-09 pro- files, as the debris were characterized by low coverage and few or no species of lichens. The debris likely originated from recent notches formed in the rockwall above the SW-08 (Fig. 10a) and NE-09 profiles (Fig. 10b). Along the other profiles (SW-01, SW-09, NE-02, NE-07), the freshest debris are found in the proximal part of the talus, presumably supplied by small-scale rockfalls. The proximal part of talus slopes generally consists of smaller debris, whereas larger debris and well-developed vegetation are found on the distal part of most of the slopes. Thus, small-scale rockfalls are unlikely to reach the base of the talus, where the accumulation is associated with large-scale rockfalls that have occurred in the past. Nevertheless, debris-size-sorting on the talus slopes, with increasing debris size toward the base, could either be attributable to rockfalls [START_REF] Kirkby | Surface stone movement and scree formation[END_REF][START_REF] Statham | A scree slope rockfall model[END_REF][START_REF] Church | Contemporary sedimentary environments on Baffin Island, NWT, Canada: debris slope accumulations[END_REF] or to talus reworking by snow avalanches [START_REF] Rapp | Recent development of mountain slopes in Kärkevagge and surroundings, northern Scandinavia[END_REF][START_REF] Luckman | Geomorphic work of snow avalanches in the canadian Rocky Mountains[END_REF][START_REF] Luckman | Debris accumulation patterns on talus slopes in Surprise Valley, Alberta[END_REF][START_REF] Jomelli | Comparing the characteristics of rockfall talus and snow avalanche landforms in an Alpine environment using a new methodological approach: massif des Ecrins, French Alps[END_REF][START_REF] Decaulne | Geomorphic evidence for present-day snowavalanche and debris-flow impact in the Icelandic Westfjords[END_REF]. Basalt debris are larger and exhibit a more massive shape than sedimentary rock debris. They also have a higher fall height, resulting in greater travel distances toward the distal part of talus slopes, as observed along the SW-06, SW-08, SW-09, NE-02 and NE-06 profiles.

Nowadays, expanding shrub cover in Tasiapik Valley [START_REF] Provencher-Nolet | Quantification des changements récents à l'écotone forêt-toundra à partir de l'analyse numérique de photographies aériennes[END_REF]Pelletier, 2015;[START_REF] Pelletier | Ecosystem changes across a gradient of permafrost degradation in subarctic Québec (Tasiapik Valley, Nunavik, Canada)[END_REF] could suggest a recent decline in debris supply, as rockfalls would not be frequent enough to limit shrub expansion on some slopes. On the NE side, a dense shrub cover (Fig. 10c) separates perched and basal talus slopes, while shrubs established themselves at the very base of the talus slope along the SW-07, SW-08 and SW-09 profiles (Fig. 10d). In both cases, shrubs have developed over highly weathered rock debris or rocky outcrops. Such a pattern of vegetation colonization highlights areas that lie beyond the reach of the most recent slope dynamics.

Slope debris redistribution

Scree slopes generally have a mean angle between 25 • and 30 • , as reported by [START_REF] Sauchyn | Particle size and shape variation on alpine debris fans, Canadian Rocky Mountains[END_REF] and [START_REF] Francou | Analysis of the segmentation in the profile of Alpine talus slopes[END_REF]. Steeper scree slopes (>30 • ) have been studied in Québec, namely in Schefferville [START_REF] Andrews | The development of scree slopes in the English Lake District and central Quebec-Labrador[END_REF], in Gaspésie [START_REF] Andrews | The development of scree slopes in the English Lake District and central Quebec-Labrador[END_REF][START_REF] Hétu | Le litage des Éboulis Stratifiés Cryonivaux en Gaspésie (Québec, Canada): Rǒle de la Sédimentation Nivéo-Éolienne et des Transits Supranivaux[END_REF][START_REF] Hétu | Effects of environmental change on scree slope development throughout the postglacial period in the Chic-Choc Mountains in the northern Gaspé Peninsula, Québec[END_REF][START_REF] Germain | Hillslope processes and related sediment fluxes on a finegrained scree slope of Eastern Canada. In: Source-to-Sink Fluxes in Undisturbed Cold Environments[END_REF]) and on the central islands of Wiyâshakimî Lake [START_REF] Decaulne | Geomorphic evidence of Holocene slope dynamics on the Canadian shield-a study from Lac à l'Eau-Claire, western Nunavik[END_REF]. Scree slopes also have a segmented profile with a steeper gradient in the proximal part and a strong basal concavity in the distal part; this slope geometry was observed on several talus slopes (Fig. 3). Given these findings, debris redistribution must be an ongoing process [START_REF] Church | Contemporary sedimentary environments on Baffin Island, NWT, Canada: debris slope accumulations[END_REF][START_REF] Francou | Analysis of the segmentation in the profile of Alpine talus slopes[END_REF]) that is likely the result of snow-avalanche rework. Their impact results in a concave inflection of the talus slopes and an increase in debris size towards the base, and can be observed on most talus slopes [START_REF] Rapp | Recent development of mountain slopes in Kärkevagge and surroundings, northern Scandinavia[END_REF][START_REF] Luckman | The geomorphic activity of snow avalanches[END_REF][START_REF] Luckman | Geomorphic work of snow avalanches in the canadian Rocky Mountains[END_REF][START_REF] Luckman | Debris accumulation patterns on talus slopes in Surprise Valley, Alberta[END_REF][START_REF] Jomelli | Comparing the characteristics of rockfall talus and snow avalanche landforms in an Alpine environment using a new methodological approach: massif des Ecrins, French Alps[END_REF][START_REF] Decaulne | Geomorphic evidence for present-day snowavalanche and debris-flow impact in the Icelandic Westfjords[END_REF]. The general concavity of talus slopes suggests that surficial debris redistribution is more important than debris supply from the rockwall. Instead, smaller debris tend to be trapped in the numerous surficial cavities that are found on the talus [START_REF] Statham | A scree slope rockfall model[END_REF][START_REF] De Blasio | Properties of talus slopes composed of flat blocks[END_REF].

Photographic monitoring of a portion of the SW side found evidence of discrete rockfall events in June and July 2018, but it has mainly highlighted numerous wet and dirty snow avalanches that occurred in the spring (April-June 2018). The latter have a greater erosive capacity due to higher friction at their base. They can also dislodge rock material from the bare rockwall, thus supplying the talus slope with new debris (Gardner, 1983a;[START_REF] Jomelli | Dépôts d'avalanches dans les Alpes franc ¸ aises: géométrie, sédimentologie et géodynamique depuis le Petit Age Glaciaire[END_REF][START_REF] Mcclung | The Avalanche Handbook[END_REF]. During the same period, the snow-covered talus enabled recently fallen debris to slide down from the apex to mid-slope (∼30 m). Sedimentary rock debris are prone to efficient sliding because of their flattened shape (the mean flatness index for sedimentary rock debris (3.04) is higher than basalt debris (2.49) on the investigated slopes), as reported by Pérez (1989) at Lassen Peak, California, and[START_REF] Hétu | Le litage des Éboulis Stratifiés Cryonivaux en Gaspésie (Québec, Canada): Rǒle de la Sédimentation Nivéo-Éolienne et des Transits Supranivaux[END_REF] in Gaspésie, Québec.

Results of reach angle measurements have shown that scattered slope debris located far downslope originated from large-scale rockfall events along most of the talus slopes, as the angle is greater than 26 • [START_REF] Corominas | Integrated landslide susceptibility analysis and hazard assessment in the principality of Andorra[END_REF]. Furthermore, low shadow angles, especially along the SW-07 and SW-09 profiles, provide evidence that their deposition is not entirely due to the rockfall itself, and that the transition from the talus to the valley floor requires an external transport agent [START_REF] Domaas | Geometrical Methods of Calculating Rockfall Range[END_REF]. The lowest measured shadow angle was 17 • in Norway [START_REF] Domaas | Geometrical Methods of Calculating Rockfall Range[END_REF], whereas it generally ranges between 22 • and 30 • [START_REF] Rapp | Recent development of mountain slopes in Kärkevagge and surroundings, northern Scandinavia[END_REF][START_REF] Govi | Photo-interpretation and mapping of the landslides triggered by the Friuli earthquake (1976)[END_REF][START_REF] Lied | Rockfall problems in Norway[END_REF][START_REF] Hungr | Engineering evaluation of fragmental rockfall hazard[END_REF][START_REF] Evans | The assessment of rockfall hazard at the base of talus slopes[END_REF]. Debris falling onto an ice-covered talus slope are unlikely to be deposited at such a great distance, considering that [START_REF] Evans | The assessment of rockfall hazard at the base of talus slopes[END_REF] reported a 24 • shadow angle for a boulder falling onto a smooth glacier. In addition, as mentioned above, the current position of the farthest lying boulders does not support the hypothesis of snow-avalancherelated transport. The abrupt slope transition (steep rockwall to valley floor) is not conducive to extreme snow-avalanche runout distances [START_REF] Bakkehoi | Calculation of snow avalanche runout distance[END_REF][START_REF] Mcclung | The Avalanche Handbook[END_REF]. From the base of the talus, the slope angle decreases to ∼4 • , which means that a boulder that is located beyond this zone could not have been transported by a snow avalanche, being located outside of the runout zone.

Therefore, the deposition of debris would most likely be related to the deglaciation/postglacial marine episode, starting at around 8200 cal. BP in the area [START_REF] Hillaire-Marcel | La déglaciation et le relèvement isostatique sur la côte est de la baie d'Hudson[END_REF][START_REF] Lavoie | Deglaciation landforms and C-14 chronology of the Lac Guillaume-Delisle area, eastern Hudson Bay: a report on field evidence[END_REF]. In particular, the presence of highly altered basalt boulders (Fig. 11) overlying littoral marine sediments would suggest that their depo-sition occurred during the regression of the Tyrrell Sea. The great distance from the surrounding slopes could then be attributable to sea-ice-related processes, such as ice rafting or ice pushing. Finally, similar processes could have distributed glacially transported rounded boulders against the talus slopes in the lower part of the NE slope, where the vast majority of the subrounded and/or gneiss debris were observed.

Rockwall erosion

Different dismantling mechanisms affect the various volcanosedimentary lithologies exposed on the cuesta frontslope. The top basalt layer shows a distinct columnar polygonal jointing, resulting in the detachment and subsequent fall of large (>5 m a-axis) monoliths (Fig. 12a). The edge of the basalt layer on the SW side reveals a sawtooth shape (Fig. 12b), also documented by [START_REF] Belzile | Les versants rocheux périglaciaires à la presqu'île des Manitounouc Kuujjarapik[END_REF] at the Manitounuk Peninsula, 100 km south of Umiujaq (55 • 42 N, 77 • 07 W). This pattern highlights the numerous monolith falls that occurred in the past. Ongoing periglacial processes such as gelifraction and frost heave caused extensive basalt dismantling (Fig. 12c). [START_REF] Michaud | Altération des substrats rocheux et rôle du soulèvement gélival dans la formation des champs de blocaille, en Hudsonie[END_REF] observed similar periglacial weathering in the basalt bedrock about 45 km south of Umiujaq (56 • 09 N, 76 • 36 W), resulting in block field development. Exhaustion of the underlying sedimentary rock layers could also cause the basalt layer to overhang the slope, eventually leading to rockfalls (Fig. 12d).

The underlying sedimentary rock layers are prone to frost shattering due to their thin-bedded sub-horizontal structure and abundant fractures exposed on the rockwall, which means that these layers are also subject to rockfalls [START_REF] Frayssines | Failure mechanisms and triggering factors in calcareous cliffs of the Subalpine Ranges (French Alps)[END_REF][START_REF] Mateos | Freeze-thaw cycles and rainfall as triggering factors of mass movements in a warm Mediterranean region: the case of the Tramuntana Range (Majorca, Spain)[END_REF][START_REF] Letortu | Le recul des falaises crayeuses hautnormandes et les inondations par la mer en Manche centrale et orientale: de la quantification de l'aléa à la caractérisation des risques induits[END_REF][START_REF] D'amato | Influence of meteorological factors on rockfall occurrence in a middle mountain limestone cliff[END_REF]. Notches have recently developed in the sedimentary rock layers, as shown on Fig. 10a-b, and have formed prominent debris cones. Subsequently, these notches will induce more discrete rockfalls from the overlying layers. Locally, rocky outcrops reveal small escarpments (∼1-5 m high) that supply debris for nearby talus through short fall height and distances, especially on the NE side.

Freeze-thaw weathering is known to be a major factor in rockfall triggering, especially on deglaciated rockwalls [START_REF] Matsuoka | Rockfall activity from an alpine cliff during thawing periods[END_REF][START_REF] Ballantyne | Paraglacial geomorphology[END_REF][START_REF] Matsuoka | Frost weathering and rockwall erosion in the southeastern Swiss Alps: long-term (1994-2006) observations[END_REF]. The air temperature data record for the Umiujaq region [START_REF] Cen | Climate station data from the Umiujaq region in Nunavik[END_REF], from September 2002 to June 2018, suggests that freeze-thaw cycles are on average more frequent during the spring (27 -March to June) than in autumn (19 -September to December), but with important interseasonal differences. For instance, only 9 cycles occurred in autumn 2007, followed by 34 cycles during the next spring. However, according to [START_REF] Hales | Climatic controls on frost cracking and implications for the evolution of bedrock landscapes[END_REF], the frequency of freeze-thaw cycles would not be a determining factor in rockfall triggering. Instead, they propose a temperature threshold (-3 • C to -8 • C) whereby thin water layers can infiltrate and cause segregated ice to form in the fractures, leading to the subsequent dismantling of the bedrock. The air temperature data from October 2002 to May 2018 (total of 5688 days) suggests that there were 642 days that had a mean temperature between -3 and -8 • C. They are more frequent in November (31.8%), April (20.9%), December (15.7%) and May (13.9%). At a depth of 100 cm in the basalt bedrock, the number of days within this threshold temperature range were more frequent in April (39.5%), December (35.2%) and January (11.9%) between 2001 and 2006 [START_REF] Allard | Borehole and near-surface ground temperatures in northeastern Canada[END_REF], which shows a potential delay from autumn to winter compared to air temperature. It should be noted that slope orientation and altitude could influence the temperature at the surface of the rockwall, and that significant differences could be expected between the SW (NE facing) and NE sides of the valley. Nevertheless, the available freeze-thaw cycle data and temperature threshold data concur with the seasonal trend (autumn/spring) for potential rockfall activity (Gardner, 1983b;[START_REF] Matsuoka | Rockfall activity from an alpine cliff during thawing periods[END_REF][START_REF] Hales | Climatic controls on frost cracking and implications for the evolution of bedrock landscapes[END_REF][START_REF] Matsuoka | Frost weathering and rockwall erosion in the southeastern Swiss Alps: long-term (1994-2006) observations[END_REF].

On a longer time scale, paraglacial adjustment is another rockfall-triggering factor that should be considered on deglaciated slopes [START_REF] Ballantyne | Paraglacial Slope Adjustment and Resedimenfation following Recent Glacier Retreat, Fåbergstølsdalen, Norway[END_REF][START_REF] Ballantyne | Paraglacial geomorphology[END_REF][START_REF] Kellerer-Pirklbauer | Paraglacial slope adjustment since the end of the Last Glacial Maximum and its long-lasting effects on secondary mass wasting processes: Hauser Kaibling, Austria[END_REF][START_REF] Cossart | An overview of the consequences of paraglacial landsliding on deglaciated mountain slopes: typology, timing and contribution to cascading fluxes. Quaternaire[END_REF][START_REF] Grämiger | Beyond debuttressing: mechanics of paraglacial rock slope damage during repeat glacial cycles[END_REF]. Deglaciation occurred at about 8200 cal. BP in the study area and induced rapid rates of isostatic uplift [START_REF] Hillaire-Marcel | La déglaciation et le relèvement isostatique sur la côte est de la baie d'Hudson[END_REF]Lajeunesse and Allard, 2003a;[START_REF] Lavoie | Deglaciation landforms and C-14 chronology of the Lac Guillaume-Delisle area, eastern Hudson Bay: a report on field evidence[END_REF]. This in turn caused an adjustment of the rock substratum and led to stressrelease on the rock slopes resulting from the ice retreat, also known as 'debuttressing' [START_REF] Ballantyne | Paraglacial geomorphology[END_REF]. The rockwall became susceptible to slope failures as a result of these changes.

Slope development results (Ho/Hi index) suggested that the SW side of the valley is at a younger stage of development, considering the height of the remaining rockwall compared to the step-like topography of the NE side. However, well-developed vegetation on the debris observed throughout the SW profiles (except for the SW-08 profiles) raises questions about the spatial distribution and frequency of the slope processes. The exhaustion model, in which slope failures occurring after ice retreat gradually decline until complete exhaustion of the rockwall [START_REF] Cruden | Exhaustion and steady state models for predicting landslide hazards in the Canadian Rocky Mountains[END_REF], would not fit in this case since local characteristics control the triggering of slope processes (i.e. periglacial processes). A bimodal/multimodal model [START_REF] Krautblatter | Glacier-and Permafrost-Related Slope Instabilities[END_REF] would thus be more representative considering that the region has undergone periods of climate fluctuations since its deglaciation (i.e. Little Ice Age). However, the occurrence of extreme meteorological events, such as warm summers [START_REF] Gruber | Permafrost thaw and destabilization of Alpine rock walls in the hot summer of 2003[END_REF] and heavy rainfalls [START_REF] Rapp | Slope erosion due to extreme rainfall in the Scandinavian mountains[END_REF][START_REF] Delonca | Statistical correlation between meteorological and rockfall databases[END_REF] should also be considered, as these events can trigger rockfalls and contribute to the slope development. Furthermore, a winter/spring period with many snow avalanches could result in increased erosion on the slope, thus favoring rockfalls, and cleaning available debris accumulated in the rockwall chutes.

Conclusion

In this study, we provided evidence of slope activity resulting in talus formation in Tasiapik Valley. The results of the topographic, granulometric, morphometric, petrographic and vegetation surveys suggest that talus slopes throughout the study area are at different stages of development, with some ancient and recent slope deposits. In the field, this results in steep and concave talus slopes with fresh debris, but also coarse deposits with an openwork texture showing alteration and a more developed vegetation cover. In addition, the opposite SW and NE sides exhibit significant differences with respect to most of the parameters surveyed, which means that their evolution occurred at different time scales.

Following the last deglaciation, paraglacial adjustment could have enabled slope failures, initiating talus slope formation. Nowadays, evidence of periglacial processes (gelifraction, frost heave) that have interacted with extensive bedrock jointing has been highlighted. As a result, features such as monoliths that have detached from the basalt layer and are on the verge of falling, or notches within the sedimentary rocks, show that rockfalls could occur at any given time. However, dense vegetation cover, primarily in the form of dense shrubs, has developed on some of the slopes and at the base of the talus slopes, meaning that the slope process runout is limited at the present time.

The present-day slope dynamic has been documented with the use of automatic cameras over a one-year period during the 2017-2018 winter season. The reworking of snow avalanches on slope deposits appeared to be a significant factor in the redistribution of debris, especially in the spring when the wet snow avalanche deposits mostly consisted of dirty snow. Discrete rockfalls have occurred during the same period, and some of the fallen debris have been transported down the snow-covered talus slope.

Fig. 1 .

 1 Fig. 1. Location of Tasiapik Valley within the Umiujaq area (A); regional geology and quaternary sediments in Tasiapik Valley (B); distribution of the talus slopes and investigated slopes (C). Sources of background images: MRNF (A, B), UMI orthomosaic (2010) (C).

Fig. 2 .

 2 Fig. 2. Location of the cameras along the SW side of Tasiapik Valley. Frame view of each camera is shown on the left. Source of background image: UMI orthomosaic (2010).

Fig. 3 .

 3 Fig. 3. Longitudinal cross-section of the investigated slopes. Source of background image: UMI orthomosaic (2010).

Fig. 4 .

 4 Fig. 4. Vegetation stages along the talus slopes.

Fig. 5 .

 5 Fig. 5. Slope development stage according to three parameters: 1) vegetation; 2) Ho/Hi index and 3) a combination of the two parameters. Source of background image: UMI orthomosaic (2010).

Fig. 6 .

 6 Fig. 6. Proportion of the different lithologies along the talus slope.

Fig. 7 .

 7 Fig. 7. Distribution of the debris size (a-axis) along the talus slopes. Close-up of the largest sampled debris along SW-09 profile (A) and NE-09 profile (B). Source of images: UMI orthomosaic (2010).

Fig. 8 .

 8 Fig. 8. Cross-section of the slopes (talus/rockwall) showing the reach and shadow angle for the farthest slope debris, and the hypothetical travel distance from the source area. Source of background image: MRFP.

Fig. 9 .

 9 Fig. 9. Snow-avalanche event exhibiting a dirty deposit, as the result of incorporated debris (A); discrete rockfall event and debris sliding onto the snow-covered talus (B).

Fig. 10 .

 10 Fig. 10. Recent notches formed in the rockwall along SW-07 (A) and NE-09 profiles (B); dense shrub cover at the base of the slope near SW-07 and SW-08 profiles (C) and mid-slope near NE-05 and NE-06 profiles (D).

Fig. 11 .

 11 Fig. 11. Cluster of basalt boulders located ∼250 m away from the slope, overlying littoral sediments.

Fig. 12 .

 12 Fig. 12. Erosion of the basalt layer at the top of the cuesta on the SW side, with tall monolith on the verge of falling (A), resulting in a sawtooth shape exposed on the rockwall (B); frost heave and extensive jointing of the basalt bedrock (C); basalt and sandstone layers overhanging quartz arenite layers on the NE side (D). Photos: Decaulne, 2017.

Table 1

 1 Topographic parameters of the investigated slopes.

	Profiles	Slope angle ( • )	Mean angle ( • )	Median angle ( • )	Inflection
	SW-01	18.4			linear
	SW-02	28.1			slightly concave
	SW-03	28.3			slightly concave
	SW-04	19.3			concave
	SW-05	24.9	25.3	26.3	concave
	SW-06	27.1			linear
	SW-07	26.3			concave
	SW-08	30.0			concave
	SW-09	25.2			complex
	NE-01	14.2			linear
	NE-02	12.3			slightly concave
	NE-03	17.5			slightly concave
	NE-04	16.1			complex
	NE-05	17.4	21.2	17.5	linear
	NE-06	31.8			linear
	NE-07	27.1			convex
	NE-08	23.2			linear
	NE-09	31.6			slightly concave

Table 2

 2 Calculation of Ho/Hi index.

	Profiles	Talus base elev. (m)	Talus apex elev. (m)	Ho value	Rockwall elev. (m)	Hi value	Ho/Hi index
	SW-01	171.3	187.6	16.3	222	50.7	0.321
	SW-02	147.3	178.3	31.0	216	68.7	0.451
	SW-03	146.3	176.8	30.5	208	31.7	0.494
	SW-04	151.4	183.2	31.8	208	56.6	0.562
	SW-05	148.0	187.9	39.9	212	64.0	0.623
	SW-06	127.4	159.2	31.8	212	84.6	0.376
	SW-07	46.5	82.0	35.5	214	167.5	0.212
	SW-08	47.6	67.2	19.6	222	174.4	0.112
	SW-09	49.6	75.9	26.3	272	222.4	0.118
	NE-01	156.1	169.9	13.8	202	45.9	0.301
	NE-02	157.2	167.3	10.1	252	94.8	0.107
	NE-03	164.6	182.8	18.2	254	89.4	0.204
	NE-04	163.0	184.3	21.3	270	107.0	0.199
	NE-05	162.1	178.1	16.0	282	119.9	0.133
	NE-06	225.6	260.8	35.2	282	54.4	0.624
	NE-07	131.7	167.5	35.8	292	160.3	0.223
	NE-08	134.6	184.0	49.4	302	167.4	0.295
	NE-09	161.2	188.8	27.6	312	150.8	0.183

Table 3

 3 Debris size and morphometric parameters along the longitudinal profiles.

	Profiles		a-axis (cm)	b-axis (cm)	c-axis (cm)	b/a	c/b	Flatness	Lengthening	Sphericity
		median	45.00	32.00	19.00	0.70	0.63	1.96	1.43	0.65
		mean	65.54	45.74	27.54	0.69	0.64	2.36	1.53	0.65
	SW-01	standard dev. maximum	58.57 450.00	44.99 350.00	26.57 180.00	0.16 0.98	0.22 1.00	1.24 10.42	0.40 2.73	0.12 0.93
		minimum	11.00	7.00	3.00	0.37	0.11	1.10	1.03	0.38
		range	439.00	343.00	177.00	0.61	0.89	9.31	1.70	0.55
		median	60.00	37.00	21.00	0.63	0.64	2.10	1.58	0.62
		mean	78.18	48.93	30.99	0.64	0.63	2.55	1.68	0.62
	SW-02	standard dev. maximum	80.30 600.00	48.04 350.00	35.48 290.00	0.17 1.15	0.21 1.04	2.44 27.50	0.51 3.64	0.11 0.87
		minimum	8.00	5.00	1.00	0.28	0.05	1.17	0.87	0.26
		range	592.00	345.00	289.00	0.88	0.99	26.33	2.77	0.61
		median	84.00	50.00	23.00	0.67	0.57	2.21	1.50	0.60
		mean	154.48	95.20	60.77	0.67	0.58	2.84	1.62	0.62
	SW-03	standard dev. maximum	262.49 1900.00	141.39 1000.00	108.40 800.00	0.18 1.00	0.24 1.11	1.88 13.75	0.49 3.32	0.13 0.91
		minimum	8.00	7.00	3.00	0.30	0.10	1.07	1.00	0.32
		range	1892.00	993.00	797.00	0.70	1.01	12.68	2.32	0.59
		median	50.00	36.00	17.00	0.72	0.58	2.21	1.38	0.61
		mean	77.14	54.89	29.70	0.72	0.56	2.69	1.48	0.64
	SW-06	standard dev. maximum	78.61 470.00	55.69 310.00	34.26 190.00	0.17 1.04	0.24 1.00	1.36 7.50	0.45 3.86	0.13 0.93
		minimum	5.00	5.00	3.00	0.26	0.15	1.06	0.96	0.33
		range	465.00	305.00	187.00	0.78	0.85	6.44	2.90	0.60
		median	23.50	14.00	6.00	0.63	0.50	3.00	1.58	0.53
		mean	30.08	17.27	8.24	0.62	0.89	4.22	1.81	0.54
	SW-07	standard dev. maximum	24.17 206.00	11.83 100.00	7.91 70.00	0.18 1.00	1.33 10.00	3.78 29.50	0.76 6.80	0.13 0.94
		minimum	6.00	2.00	1.00	0.15	0.04	1.10	1.00	0.21
		range	200.00	98.00	69.00	0.85	9.96	28.40	5.80	0.73
		median	36.00	20.00	11.00	0.59	0.56	2.57	1.69	0.56
		mean	48.16	27.31	14.29	0.60	0.57	2.96	1.86	0.57
	SW-08	standard dev. maximum	42.23 280.00	26.74 223.00	12.83 100.00	0.18 1.00	0.22 1.00	1.53 11.17	0.80 10.42	0.12 0.92
		minimum	4.00	3.00	1.00	0.10	0.09	1.07	1.00	0.19
		range	276.00	220.00	99.00	0.90	0.91	10.10	9.42	0.73
	measured clasts. The largest clasts range in size between 350 and
	950 cm, while the mean flatness index values range from 2.29 to
	3.58. Lengthening and sphericity indices show similar values as the
	profiles on the SW side, ranging from 1.40 to 1.82, and from 0.57
	to 0.69. By comparing the calculated index values

Acknowledgements

Funding for this project was provided by the Natural Sciences and Engineering Research Council (NSERC), LabEx DRIIHM and OHMi NUNAVIK-TUKISIG, and IPEV program DeSiGN. The authors want to thank Félix Faucher, Julien Lebrun and Thorsteinn Saemundsson for their valuable help in the field and the Centre d'études nordiques (CEN) for its logistical support.


R Razzaghi

F Colas

X Guillaud

M Paolone

F Rachidi

Hardware-in-the-Loop Validation of an FPGA-Based Real-Time Simulator for Power Electronics Applications

Keywords: Hardware-in-the-loop, real-time simulation, field programmable gate array, modified nodal analysis, fixed admittance matrix nodal method

This paper presents the hardware-in-the-loop (HIL) validation of a proposed FPGA-based real-time simulator for power electronics applications. The proposed FPGA-based realtime simulation platform integrates the Modified Nodal Analysis (MNA) method, Fixed Admittance Matrix Nodal Method (FAMNM) and an optimization technique to assess the optimal value of the switches conductance in order to minimize the relevant errors. Moreover, the proposed platform includes an automatic procedure to translate the netlist user-defined circuit schemes to the relevant equations to be solved in the FPGA. Then, the paper illustrates the validation of the proposed simulator in two steps. First the validation is presented by comparing the FPGA-based simulation results with the offline ones performed by EMTP-RV. Then, further validation is presented by means of a dedicated HIL experimental setup composed of a controller connected to an actual two-level, three-phase inverter and its corresponding FPGA real-time model.

I. INTRODUCTION

Real-time simulation is a way to couple replica models of a given hardware, or system, with real-scale monitoring and control devices/systems. Such simulations are referred as hardware-in-the-loop (HIL) ones and allow performing different operational or control experimental tests which might not be possible to be conducted on the real hardware/system (e.g., [START_REF] Lu | A low-cost real-time hardware-in-the-loop testing approach of power electronics controls[END_REF], [START_REF] Ould-Bachir | Effective floating-point calculation engines intended for the FPGAbased HIL simulation[END_REF]). HIL simulations, based on the exchange of analog/digital signals between the hardware under test (HUT) and the simulator, are categorized as controller HIL (CHIL) or power HIL (PHIL) (e.g., [START_REF] Ren | Improve the Stability and the Accuracy of Power Hardware-in-the-Loop Simulation by Selecting Appropriate Interface Algorithms[END_REF]). Concerning the CHIL, the HUT is a controller and the exchanged signals usually have low amplitudes. PHIL simulations refer to a kind of HIL tests where the HUT is a power device and, typically, proper amplifiers are required to adapt the exchanged physical signals.

For industrial applications, there are two main types of hardware used to develop a real-time simulator for the HIL tests: (i) CPU-based simulators, (ii) FPGA-based ones. In general, CPU based real-time simulators represent a better R. Razzaghi and M. Paolone are with the Distributed Electrical Systems Laboratory, F. Rachidi is with the Electromagnetic Compatibility Laboratory. (EMC). Both laboratories are with the École Polytechnique Fédérale de Lausanne (EPFL), Switzerland, (e-mail: reza.razzaghi@epfl.ch, phone: +41789075556, mario.paolone@epfl.ch, farhad.rachidi@epfl.ch). X. Guillaud and F. Colas are members of the L2EP laboratory. X. Guillaud is with Ecole Centrale de Lille, France (e-mail: xavier.guillaud@ec-lille.fr). F. option to simulate bulk power networks since they can achieve acceptable simulation time steps (e.g., in the order of few tens of microseconds) and represent relatively complex systems. Additionally, existing CPU-based real-time simulators are typically linked to well-established programming environments (e.g., MATLAB SimPowerSystems (SPS)) that allow a more straightforward way to model components and run the simulation. However, the achievable integration time steps of CPU-based real-time simulators have a lower bound associated with the partial sequential operations that the CPU architectures need to deploy. As a consequence, the relatively large simulation time steps required by these simulators do not allow to model high frequency phenomena such as electromagnetic transients in power converters. With particular reference to this last item, as indicated in [START_REF] Gole | Guidelines for modeling power electronics in electric power engineering applications[END_REF], the simulation time-step should be at least 20 times smaller than the switching frequency. Therefore, further techniques (e.g., interpolation ones) are required to be employed to qualify CPU-based simulators for high PWM power electronics [START_REF] Dufour | Review of state-of-the-art solver solutions for HIL simulation of power systems, power electronic and motor drives[END_REF].

During the past years, the size and computational power of FPGAs has been dramatically increased. As a consequence, FPGA-based real-time simulation has emerged as a leading trend for the Electromagnetic Transient (EMT) simulations of power systems and HIL simulations (e.g., [START_REF] Matar | Massively Parallel Implementation of AC Machine Models for FPGA-Based Real-Time Simulation of Electromagnetic Transients[END_REF], [START_REF] Dufour | General-purpose reconfigurable low-latency electric circuit and motor drive solver on FPGA[END_REF]).

Concerning the HIL simulation of power electronics applications, FPGA-based real-time simulators provide several advantages over CPU-based ones (e.g., [START_REF] Ould-Bachir | Effective floating-point calculation engines intended for the FPGAbased HIL simulation[END_REF], [START_REF] Dufour | Review of state-of-the-art solver solutions for HIL simulation of power systems, power electronic and motor drives[END_REF]). In particular, the parallel processing hardwired in FPGAs enables the implementation of specific methodologies that dramatically reduce the sequencing of the operations taking place in CPUs. FPGA-based real-time simulators provide lower sampling rate, higher frequency bandwidth and lower I/O latency [START_REF] Dufour | Review of state-of-the-art solver solutions for HIL simulation of power systems, power electronic and motor drives[END_REF].

However, the FPGA-based real-time simulations are characterized by some limitations. In particular, the model development requires, in general, the knowledge of the Hardware Description Language (HDL) which limits the implementation of the complex models.

Moreover, the matrix manipulation operations are limited in the FPGAs and, as a consequence, the simulation of the switching devices such as power electronics requires particular treatments. In this respect, the most straightforward method to represent topology-variable circuits in FPGA real-time simulators is the so-called Fixed Admittance Matrix Nodal Method (FAMNM) [START_REF] Dufour | General-purpose reconfigurable low-latency electric circuit and motor drive solver on FPGA[END_REF]. This method, irrespective of the number of the switches and their states, allows for obtaining a fixed nodal admittance matrix during switching transitions. However, it introduces artificial oscillations and errors in the simulation results (e.g., [START_REF] Razzaghi | A method for the assessment of the optimal parameter of discrete-time switch model[END_REF]).

Recently, several studies have been performed in the literature on the applicability of the FPGA-based real-time simulators for the HIL simulation of power electronics applications (e.g., [START_REF] Ould-Bachir | Effective floating-point calculation engines intended for the FPGAbased HIL simulation[END_REF], [START_REF] Dufour | General-purpose reconfigurable low-latency electric circuit and motor drive solver on FPGA[END_REF], [START_REF] Matar | FPGA Implementation of the Power Electronic Converter Model for Real-Time Simulation of Electromagnetic Transients[END_REF], [START_REF] Dagbagi | FPGA implementation of Power Electronic Converter real-time model[END_REF], [START_REF] Razzaghi | A general purpose FPGAbased real-time simulator for power systems applications[END_REF]). These studies, mainly, are based on the use of FMANM approach which enables simulation of power electronics within very low time steps. However, they do not take into account the tuning of the discrete-time switch conductance value and its effect on the simulation results accuracy. Moreover, the obtained FPGAbased real-time simulation results are validated by comparing them with the offline simulation ones (e.g., SPS or EMTP-RV).

Within this context, this paper briefly illustrates a method to develop FPGA-based real-time simulation for power electronics applications that integrates the Modified Nodal Analysis (MNA) method, FAMNM and an optimization technique proposed in [START_REF] Razzaghi | A method for the assessment of the optimal parameter of discrete-time switch model[END_REF] to find the optimal value of the switches conductance in order to minimize the relevant errors. The proposed method includes an automatic procedure to translate the netlist user-defined circuit schemes to the relevant equations to be solved in the FPGA. Then, the paper mainly focuses on illustrating the validation of the proposed simulator by means of a dedicated HIL experimental setup composed of a controller connected to an actual two-levels, three-phase inverter and its corresponding FPGA real-time model.

The structure of this paper is as following: Section II provides a brief overview of the EMT simulation. Section III describes the proposed real-time simulation platform. Section IV illustrates the experimental HIL setup, the FPGA model of the two-level three-phase inverter, its preliminary validation by comparing its results with the offline simulation ones and comparison with an actual three-phase inverter. Section VI concludes the paper with final remarks.

II. OVERVIEW OF EMT SIMULATIONS

A. Circuit solvers and numerical integration methods

In general, two main types of solution methods are used in power systems and power electronics electromagnetic simulations: (i) nodal and (ii) state-space ones [START_REF] Mahseredjian | Computation of power system transients: overview and challenges[END_REF]. In this paper we have adopted the first one since it allows straightforward formulation of the power electronics systems equations and, in particular, it enables the FAMNM approach. MNA is represented by the general equation of (1) [START_REF] Chua | Computer-Aided Analysis of Electronic Circuits: Algorithms & Computational Techniques[END_REF]:

[ ][ ] [ ] n n n A x b = (1)
where matrix [A n ] is formed by the discrete representation of the network elements; [x n ] is the vector of unknowns including the network's node voltages and branch currents; and [b n ] is a vector of the independent sources and current history terms related to the network components. For EMT simulation applications, trapezoidal and backward-Euler methods are the most popular numerical integration methods [START_REF] Mahseredjian | Computation of power system transients: overview and challenges[END_REF]. For the case of switching devices, it is preferred to use the latter one since backward-Euler rule gives better damping to numerical oscillations introduced by switches [START_REF] Pejovic | Method for fast time-domain simulation of networks with switches[END_REF].

B. Discrete models of simple network components

1) Lumped elements (L,C) The most common approach for discrete-time representation of the network elements is the one proposed in [START_REF] Dommel | Digital Computer Solution of Electromagnetic Transients in Single-and Multiphase Networks[END_REF] where the circuit elements are converted into their Norton equivalent. In particular, the lumped elements (R, L, C) connected between nodes k and m are described by [START_REF] Chua | Computer-Aided Analysis of Electronic Circuits: Algorithms & Computational Techniques[END_REF], [START_REF] Dommel | Digital Computer Solution of Electromagnetic Transients in Single-and Multiphase Networks[END_REF]:

( ) ( ) ( ) ( ) ( ) eq k m km hist G v t v t i t I t t - = + -Δ (2)
where eq G is the equivalent conductance, and

( )

hist I t t -Δ is the
history current source associated with the time-discretized element. The values for the equivalent conductance and the history current are determined by the element type (i.e., R, L, C) together with the adopted numerical integration method [START_REF] Dommel | Digital Computer Solution of Electromagnetic Transients in Single-and Multiphase Networks[END_REF].

2) Switches Accurate and efficient switch modeling is a challenging issue for EMT simulators, especially when real-time constraints need to be achieved. In general, detailed switch models are too much sophisticated and not suitable for real-time applications. Therefore, behavioral switch models have been proposed for the EMT real-time applications [START_REF] Jin | Behavior-mode simulation of power electronic circuits[END_REF]. Among them, the simplest ones are the ideal switch model or the so-called two-valued resistor model where two resistors, characterized by large differences of their resistances, are associated with each state of the switch (R off , R on ). However, as well described by the literature on the subject, for both models the system's admittance matrix needs to be updated and re-factorized after each switching change generating major issues to satisfy the FPGA computational time constraints.

On the contrary, the use of the discrete-time switch model allows defining the so-called fixed nodal admittance matrix method (FAMNM). In this case, the switch is represented by a relatively small inductance when its state is 'closed' and by a relatively small capacitance when its state is 'open' (e.g., [START_REF] Pejovic | Method for fast time-domain simulation of networks with switches[END_REF], [START_REF] Hui | A discrete approach to the modeling of power electronic switching networks[END_REF]). As a consequence, in view of (2), the switch is replaced by an equivalent conductance (G s ) in parallel with a controlled current source.

The main drawback of such representation is that it introduces artificial parameters in the circuit and, consequently, oscillations and errors to the results [START_REF] Razzaghi | A method for the assessment of the optimal parameter of discrete-time switch model[END_REF], [START_REF] Mu | Improved ADC Model of Voltage-Source Converters in DC Grids[END_REF]. Therefore, an optimal tuning of the switch conductance value is an important issue to achieve accurate results. An efficient method for the optimal selection of this parameter has been proposed in [START_REF] Razzaghi | A method for the assessment of the optimal parameter of discrete-time switch model[END_REF] and it is the method adopted in this paper to properly select this parameter (see [START_REF] Razzaghi | A method for the assessment of the optimal parameter of discrete-time switch model[END_REF] for further details).

III. THE PROPOSED REAL-TIME SIMULATION PLATFORM

The overall structure of the proposed and developed real-time simulator is schematically represented in Fig. 1. In what follows the various blocks appearing in this figure are explained.

With reference to the adopted hardware platform, the proposed FPGA-based real-time simulator is based on the National Instruments compactRIO-9033, an industrial reconfigurable real-time embedded hardware platform combining an Intel Atom dual-core processor, a Xilinx Kintex-7 FPGA, and reconfigurable I/O modules. This embedded system is based on NI Linux Real-Time OS and is programmed by using the NI LabVIEW-FPGA environment. The reason to choose this hardware platform is that it provides reconfigurable platform including the CPU and the FPGA as well as reconfigurable I/O modules.

A. Circuit pre-processing

In the proposed real-time simulation platform, the EMTP-RV simulation environment is used as a Graphical User Interface (GUI) to define the circuit under study and its parameters. Then, the designed circuit is analyzed by this software to generate the so-called netlist file. It contains all the information about the types of the circuit components, their values, and their interconnections.

The netlist file is then used by netlist pre-processing unit to extract the relevant information to be passed to the FPGA solver. Since this is an offline process, it is done by the CPU of the real-time hardware platform.

According to the type of the element indicated in the netlist file, the algorithm extracts the relevant information (e.g., topological connections, values, etc.). The structure of the data pre-processing unit is shown in Fig. 1.

The output arrays of the netlist pre-processing unit are used by NAM builder block to form the nodal admittance matrix (NAM) in the simulator. This matrix is inverted in the CPU level based on the floating point numerical representation and the double-precision. The inverted matrix is transformed to the fixed-point numerical representation form by proper selection of the fixed point in order to provide good accuracy. Then, the matrix is stored in the memory blocks in order to be transferred to the FPGA solver.

It is worth noting that, the optimization problem to find the optimal switches conductance values is performed in the CPU of the real-time hardware platform by the G s optimization unit. Then, the calculated values are used to build the nodal admittance matrix and, also, transferred to the FPGA to be used in the switches RHS computations.

Moreover, additional data are calculated in the other parameters unit and transferred to the FPGA. These data include desired simulation time step, independent voltage/current sources information, the number of elements in the network per element type, and eventual controller variables (in our specific case, the converter PWM setup variables).

B. FPGA circuit solver

In order to take advantage of the parallel processing capability of the FPGAs, the adopted solver architecture is based on several parallel sub-tasks. In particular, in order to achieve very low simulation time steps, the two main steps of solving the MNA equation ( 1) are decoupled to several parallel and independent tasks. Among them, the RHS vector update for different elements is done independently and in parallel. Namely, dedicated RHS computation units are considered for inductors, capacitors, and switches.

1) RHS_L i , RHS_C i

For the case of inductors and capacitors, the RHS variables are function of corresponding node voltages and branch currents in the previous time step (the RHS element for the resistor is zero). Therefore, the required values to compute RHS elements are stored in the FPGA memory to be accessed in the next iteration. Then, the stored variables are used to update the RHS elements. It is worth observing that, in the CPU-based real-time simulators, the update of the RHS elements is done sequentially for different types of elements. However, thanks to the inherent parallel processing capability of the FPGA, these tasks are done in parallel. In particular, for inductors and capacitors, a dedicated computational unit has been coded. The RHS elements for the inductors and capacitors can be expressed x is the inductor (or capacitor) state variable in the previous time step. For the case of inductors,

, 1 i i L C γ = , i n n i L
x i = , and for the capacitors

, i i i L C C t γ = -Δ , i n n i C x v = .
This equation is solved independently for every capacitors and inductors to realize the highest level of parallelism.

2) RHS_Sw i

The RHS elements of the switches are calculated using another dedicated sub-module. In particular, after calculating the optimal conductance values in the offline pre-processing, these values are transferred to the FPGA to be used in this submodule. Then, according to the switches states and by accessing to their voltages and currents, the RHS elements are calculated as [START_REF] Pejovic | Method for fast time-domain simulation of networks with switches[END_REF]: G is the i th switch optimal conductance value, 1 s n+ is the switch current state, and sw N is the number of switches. Similar to the inductors and capacitors, the RHS elements for different switches are calculated independently.

1 1 1 s s 1 ; 1,..., 0 i i i i n s n sw n s s n sw n i J G v i N + + + - = = ⎧ ⎪ = ⎨ = ⎪ ⎩ (4)
In order to determine the switch state, Switch states block is considered. The switch state is determined by its type (e.g., diode, IGBT-diode pair, etc.) [START_REF] Pejovic | Method for fast time-domain simulation of networks with switches[END_REF], and the switches commands can be determined by the digital input modules (Gate signals block in Fig. 1) or user-defined logics.

For the case of an IGBT in parallel with an anti-parallel diode, the switch current state is determined based on (5):

( ) ( )

1 1 0 0 n n n n n n s c s i s v + + = + ≤ + < (5)
where andn s is the switch previous state. Concerning the matrix to vector multiplication, thanks to FAMNM switch representation, the NAM is constant during the simulation. Thus, it is computed once in the pre-processing unit. In principle, the matrix-vector multiplication process consists of two loops where the outer loop is associated with the number of the matrix rows and the inner loop corresponds to the number of elements within each row (i.e., number of columns).

1 n s + is the switch current state, 1 n c + is the switch gate command,
In order to accelerate the multiplication, different levels of parallelism and techniques can be applied. In particular, the multiplication is done by splitting the matrix into individual rows and doing the dot-product and accumulation for each row, individually. Then, within each dot product operation, the multiplication is done in parallel. To this end, NI LabVIEW FPGA IP Builder tool is used to optimize the multiplication algorithm based on the requested latency and the throughput and by considering the available FPGA resources [START_REF]Using NI LabVIEW FPGA IP Builder to Optimize and Port VIs for Use on FPGAS[END_REF].

It is worth observing that, the FPGA-based calculations are based on fixed point numerical representation. In general, floating point offers higher precision for the numerical representation compared to the fixed point one. However, fixed-point representation is more efficient from the hardware resources usage and performance point of views. By carefully selecting the fixed-point representation, good accuracy values can be achieved.

Apart from the circuit solver engine, additional logics concerning the PWM controller are implemented in FPGA in order to provide higher precision for the high frequency PWM signals. The internal PWM controller logic can be used to verify the performance of the simulator without the need of using an external one.

IV. EXPERIMENTAL HIL SETUP AND FPGA MODEL OF THE SYSTEM UNDER STUDY

A. Description of the HIL setup

In this paper, to validate the performance of the developed FPGA-based real-time simulator, a two-level three-phase inverter that is connected to an inductive filter (10 mH) and a resistive load (20 Ω) is considered. The global setup is depicted in Fig. XX. It shows three main parts: the system under study which can be a FPGA based real time model or a real inverter, the DS 1104 controller board and the HMI (Human Machine Interface). A classical dq synchronous frame current controller [START_REF] Kazmierkowski | Current control techniques for three-phase voltage-source PWM converters: a survey[END_REF] has been used and implemented in the controller board. Then, the pre-processing unit analyzes the generated netlist file and builds the corresponding NAM as equation ( 7). This matrix is inverted and converted to the fixed point representation based on 40 bits for the word length and 19 bits for the integer part. It is worth observing that the model provided by [START_REF] Dufour | Review of state-of-the-art solver solutions for HIL simulation of power systems, power electronic and motor drives[END_REF] does not explicitly appear in [START_REF] Dufour | General-purpose reconfigurable low-latency electric circuit and motor drive solver on FPGA[END_REF] since it is a logic determining the status of this aggregated IGBT+diode device.

G G G G G G G G G G G G G G ⎡ - - - - - - - - - - - - - ⎣ 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 1 1 2 1 3 1 1 1 2 1 3 1 4 1 5 1 6 1 n n n n n n n n n L n L n L n S n S n S n S n S n S n dc V V V V V V V V
+ + + + ⎤ ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎢ ⎥ ⎦ 1 1 1 1 1 1 2 1 3 1 4 1 5 1 6 _ 1 _ 2
⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ = ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎣ ⎦ (7)
The optimization process to find the optimal values for the switch conductances is performed in the pre-processing phase. By considering the switching modes where each switch conducts for 180 degrees of a cycle, there are eight possible switching permutations as: (S1,S2,S6), (S1,S2,S3), (S2,S3,S4), (S3,S4,S5), (S4,S5,S6), (S1,S5,S6), (S1,S3,S5), and (S2,S4,S6). Therefore, we obtain eight ANAM corresponding to each status of the ideal switches. Among them, six switching patterns generate non-zero voltage across the load and two of them (the upper or lower switches are conducting) generate zero voltage across the load.

According to the method presented in [START_REF] Razzaghi | A method for the assessment of the optimal parameter of discrete-time switch model[END_REF], the first step to calculate the optimal value for the switch conductance is to find the sets of eigenvalues corresponding to the possible switching permutations and also, the ones of the FAMNM. It is worth noting that, since the load and filter parameters are identical for all the phases, the eigenvalues for the two sets of patterns are equal for each set. Therefore, in view of the symmetrical nature of the circuit, one identical conductance value can be assigned to the six switches. By applying the optimization method presented in [START_REF] Razzaghi | A method for the assessment of the optimal parameter of discrete-time switch model[END_REF], the objective function exhibits the optimal value of G s =0.51 (see [START_REF] Razzaghi | A method for the assessment of the optimal parameter of discrete-time switch model[END_REF] for further details about the objective function definition). This value is used to build the NAM and also update the switches RHS elements.

In the first step, the performance of the proposed FPGAbased real-time simulator is validated by comparing its results with offline simulations carried out in EMTP-RV [REF?]. Fig. 4 illustrates the three-phase load currents obtained by the FPGA-based real-time simulator, and by EMTP-RV off-line simulation environment, respectively. In this figure, the PWM carrier frequency is 1 kHz. Fig. 5 shows the error between the load currents of the benchmark model and the ones of the FPGA simulator. The calculated error is normalized based on the current peak value and shown in perunit. It can be observed that the FPGA-based results are characterized by a maximum error of 0.0015 pu with respect to the benchmark simulation. The reasons behind this small error are two: (i) the truncation realized by the fixed-point simulation calculations, and (ii) the discrete-time switch model which, as known, introduces approximations in the discrete-time switch model. Concerning the achieved integration time step, the simulation needs 6 FPGA ticks per time step. Consequently, by considering the 40 MHz FPGA clock, it results into an integration time step of 150 ns. Therefore, the availability of a faster FPGA clock will directly enable to further reduce the integration time step. The total utilized FPGA resources (based on the fixed point calculations) are: slice registers: 10.1%, slice LUTs: 28%, block RAMs: 1.5%, DSP48s: 85.3%

C. Validation by means of HIL simulation test

In the previous section, the validation of the proposed FPGA-based real-time simulator is presented by comparing its obtained results with the offline simulation ones which exhibits excellent simulation accuracy together with very low simulation time-step. In this section, a further validation test is presented by making reference to a HIL simulation test performed by the proposed FPGA-based real-time simulator. To this end, first, the experimental setup explained in the previous section is used to perform the HIL test with the proposed real-time simulation platform. The external PWM controller (the PWM frequency is 1 kHz) is coupled with the simulator by using digital input modules. In particular, the switches gate signals are determined by NI-9401 which is a high speed digital I/O module and the gate signals loop and the simulation one are synchronized. To export the simulation generated signals, NI-9263 module is used which is an analog output one. Since, the maximum sampling rate of this module is 100 kHz, the generated load current signals are downsampled by this frequency to be monitored in the oscilloscope.

Then, the same controller is coupled with a physical inverter which is connected to the physical inductive filter and resistive load with the same value of the HIL simulation. The controller type and parameters are identical to the ones of the HIL test. The load currents are measured by using the current sensors described before and are observed by the oscilloscope. Fig. 6 shows the comparison of the HIL simulation results and the measured ones for the three phase load currents and Fig. 7 shows the relative errors between the HIL simulation results and the measured ones. It is observed that the HIL simulation results are with agreement of the measurement ones with obvious higher errors compared to the offline simulations. The reasons for this error are: (i) the presence of noise in the measurements, (ii) the limited current sensors bandwidth, (iii) the error between the simulated physical ones converter filter and load, and (iv) non-linear behavior of the switches in the real inverter compared to the linear switch model in the HIL simulation. However, is spite of the above-listed source of errors, the comparison appears satisfactory providing a good experimental validation of the proposed FPGA simulation platform.

V. CONCLUSION

The paper presented the HIL validation of the proposed FPGA-based real-time simulator for the power electronics applications. The proposed real-time simulator was implemented in an industrial real-time embedded system (National Instruments CompactRio real-time platforms) and has the following features: (i) it makes use of the Modified Nodal Analysis (MNA) method, (ii) it integrates the Fixed Admittance Matrix Nodal Method (FAMNM) together with the optimal selection of the switch conductance parameter, (iii) it enables the possibility to accurately reproduce electromagnetic switching transients taking place in power electronic switching devices, and (iv) it enables to reach extremely low integration time steps (in the order of hundreds of ns) and avoids the need of redesigning and recompiling the FPGA code.

The performance of the proposed real-time simulation platform has been validated, in the first approach, by simulation of a two-level three-phase inverter and comparing the obtained real-time simulation results with the offline benchmark performed in EMTP-RV simulation environment. Then, further validation has been done by HIL simulation test performed by the proposed FPGA-based real-time simulator and comparing the obtained results with the measured ones using the real inverter connected to the physical filter and load.

The comparison of the obtained results with respect to offline benchmark ones and the measured physical ones showed a good agreement together with high computation efficiency.

Fig. 1 .

 1 Fig. 1. The FPGA solver architecture together with the real-time processer tasks.

Fig. 2 .

 2 Fig. 2. Architecture of the netlist pre-processing unit.

FPGA

  element for the i th inductor or capacitor in the current time step, corresponding to the i th inductor (or capacitor) and n i

Jv

  + is the RHS variable for the i th switch, is the i th switch voltage, i s

Fig

  Fig. XX -global setup presentation A picture of this setup is depicted in Fig. XXX. It shows the

Fig. 3 .

 3 Fig. 3. Schematic representation of the two-level three-phase inverter.

Fig. 4 .

 4 Fig. 4. Comparison of the FPGA-based real-time solver results with the corresponding EMTP-RV ones (three-phase load currents).

Fig. 5 .

 5 Fig. 5. Relative error of the load currents in pu. (reference values of Fig. 4).

Fig. 6 .

 6 Fig. 6. Comparison of the FPGA-based HIL test results with the measured ones (three-phase load currents).

Fig. 7 .

 7 Fig. 7. Relative error of the load currents in pu. (reference values of Fig. 6).

levels, three phase inverter System under study: FPGA based real time model or Real device