Xavier Bodin
email: xavier.bodin@univ-smb.fr

Emmanuel Thibert
email: emmanuel.thibert@irstea.fr

Olivier Sanchez
email: oliviersanchez64@gmail.com

Antoine Rabatel
email: antoine.rabatel@univ-grenoble-alpes.fr

Stéphane Jaillet
email: stephane.jaillet@univ-smb.fr

remote sensing

Keywords: rock glacier deformation, very high-resolution DEM, surface displacement, image correlation, TLS, ALS, SFM-MVS

Rock glaciers result from the long-term creeping of ice-rich permafrost along mountain slopes. Under warming conditions, deformation is expected to increase, and potential destabilization of those landforms may lead to hazardous phenomena. Monitoring the kinematics of rock glaciers at fine spatial resolution is required to better understand at which rate, where and how they deform. We present here the results of several years of in situ surveys carried out between 2005 and 2015 on the Laurichard rock glacier, an active rock glacier located in the French Alps. Repeated terrestrial laser-scanning (TLS) together with aerial laser-scanning (ALS) and structure-from-motion-multi-view-stereophotogrammetry (SFM-MVS) were used to accurately quantify surface displacement of the Laurichard rock glacier at interannual and pluri-annual scales. Six very high-resolution digital elevation models (DEMs, pixel size <50 cm) of the rock glacier surface were generated, and their respective quality was assessed. The relative horizontal position accuracy (XY) of the individual DEMs is in general less than 2 cm with a co-registration error on stable areas ranging from 20-50 cm. The vertical accuracy is around 20 cm. The direction and amplitude of surface displacements computed between DEMs are very consistent with independent geodetic field measurements (e.g., DGPS). Using these datasets, local patterns of the Laurichard rock glacier kinematics were quantified, pointing out specific internal (rheological) and external (bed topography) controls. The evolution of the surface velocity shows few changes on the rock glacier's snout for the first years of the observed period, followed by a major acceleration between 2012 and 2015 affecting the upper part of the tongue and the snout.

Introduction

Because of their abundance in many poorly-glacierized mountain regions of the world [START_REF] Barsch | Indicators for the Present and Former Geoecology in High Mountain Environments[END_REF], rock glaciers are key features to understand the response of the high altitude cryosphere to climate change [START_REF] Haeberli | Permafrost creep and rock glacier dynamics[END_REF]. Rock glaciers typically flow downslope at relatively low and steady rates, exhibiting surface velocity typically ranging from a few cm to 2 m per year [START_REF] Barsch | Indicators for the Present and Former Geoecology in High Mountain Environments[END_REF]. In the Alps, synchronous inter-annual fluctuations of surface velocities have been observed on many rock glaciers, suggesting a control of the ground thermal state on the deformation rate of the permafrost [START_REF] Delaloye | Recent interannual variations of rock glacier creep in the European Alps[END_REF][START_REF] Kellerer-Pirklbauer | Interannual variability of rock glacier surface velocities and its relationship to climatic conditions on a decadal scale: Some insights from the European Alps[END_REF]. Some cases of rock glaciers' destabilization were also observed in various mountain chains, such as in the Alps [START_REF] Roer | Observations and considerations on destabilizing active rock glaciers in the European Alps[END_REF][START_REF] Bodin | Collapse of the Bérard Rock Glacier (Southern French Alps)[END_REF] or in the Andes [START_REF] Bodin | Recent collapse of rock glaciers: two study cases in the Alps and in the Andes[END_REF], becoming a potential source of hazards [START_REF] Krautblatter | Glacier-and permafrost-related slope instabilities[END_REF]. These temporal changes in rock glaciers' deformation patterns pointed out various mechanical processes involved in their behavior [START_REF] Kenner | Investigating the dynamics of a rock glacier using terrestrial laser scanning, time-lapse photography, in situ GPS measurements and satellite SAR inter-ferometry: Ritigraben rock glacier[END_REF]. In this context, developing a sound knowledge of the spatial and temporal variabilities of the surface velocity of rock glaciers, as a proxy of their deformation and an indicator of their climatic controls, is needed. It is indeed a prerequisite for future efforts aiming at analyzing their physico-mechanical behavior and the driving factors, among which the climate probably plays a substantial role [START_REF] Delaloye | Recent interannual variations of rock glacier creep in the European Alps[END_REF][START_REF] Kääb | Climate change impacts on mountain glaciers and permafrost[END_REF][START_REF] Springman | Multidisciplinary investigations on three rock glaciers in the Swiss Alps: Legacies and future perspectives[END_REF][START_REF] Kenner | Estimating the Effect of Different Influencing Factors on Rock Glacier Development in Two Regions in the Swiss Alps[END_REF].

Because investigating the internal deformation of rock glaciers is very costly and limited to a single one-dimensional vertical profile along boreholes, the surface displacement is the most easily "accessible" variable for characterizing rock glacier dynamics and deformation patterns. The kinematics of the topographical surface integrates the overall vertical profile of deformation, including creeping and shearing, and their changes in space and in time are valuable information for understanding the factors that rule the deformation. Surface kinematics of rock glaciers are generally measured and monitored either remotely, with terrestrial, aerial or spatial vectors (see, e.g., [START_REF] Kääb | Analysing the creep of mountain permafrost using high precision aerial photogrammetry: 25 years of monitoring Gruben rock glacier[END_REF][START_REF] Arenson | Detection and Analysis of Ground Deformation in Permafrost Environments: Ground Deformation in Permafrost Environments: Detection and Analysis[END_REF]), or with in situ measurements, by geodetic methods (total station, differential GPS; see, e.g., [START_REF] Lambiel | Contribution of real-time kinematic GPS in the study of creeping mountain permafrost: examples from the Western Swiss Alps[END_REF]). The first set of approaches has the advantage of a wide spatial coverage, using imagery with a pixel resolution of a few decimeters to a few meters (ground sample distance (GSD)). It is nevertheless limited by the temporal availability of datasets, typically several years to decades. In the second set of methods, the surface topography can be acquired more frequently (as often as an operator can go out in the field to carry out surveys), but it results in poorly representative data in spatial terms, because it consists of punctual information (generally a marked block at the surface that can be tracked from year-to-year). Other methods for continuous measurements and monitoring of ground surface changes, like continuous L1-GPS [START_REF] Wirz | Spatio-temporal measurements and analysis of snow depth in a rock face[END_REF] or fixed camera [START_REF] Kenner | Estimating the Effect of Different Influencing Factors on Rock Glacier Development in Two Regions in the Swiss Alps[END_REF], are allowed by very recent technological developments, but remain up to now costly and not so easy to implement.

Terrestrial and airborne laser scanning (TLS, ALS) and terrestrial photogrammetry (or structurefrom-motion/multi-views stereophotogrammetry (SFM-MVS)) constitute intermediary approaches. Because of their relative simplicity in handling, these approaches easily allow one to generate (very) high-resolution digital elevation models (DEM). In addition, their repetitivity can be annual or even at smaller time intervals, mostly depending on the terrain characteristics and the range of surface velocities. So far, these methods have rarely been used for rock glacier studies [START_REF] Piermattei | Suitability of ground-based SfM-MVS for monitoring glacial and periglacial processes[END_REF][START_REF] Dall'asta | Unmanned Aerial Systems and DSM matching for rock glacier monitoring[END_REF][START_REF] Avian | First Results on Monitorin Glacier Dynamics with the Aid of Terrestrial Laser Scanning on Pasterze Glacier (Hohe Tauern, Austria)[END_REF][START_REF] Bodin | High-Resolution DEM Extraction from Terrestrial LIDAR Topometry and Surface Kinematics of the Creeping Alpine Permafrost: The Laurichard Rock Glacier Case Study (Southern French Alps)[END_REF], and no investigations have yet addressed the potential of those methods for long-term monitoring of the surface kinematics.

Aims and Objectives

In order to investigate the multi-temporal kinematics of rock glaciers and, ultimately, to better understand their climatic controls, we took advantage of an easily-accessible site with one of the longest time-series of ground-truth data (since the early 1980s; see [START_REF] Francou | 10 year surficial velocities on a rock glacier (Laurichard, French Alps)[END_REF]). On this active rock glacier located in the French Alps, we carried out several TLS, ALS and SFM-MVS surveys. The objectives of the present study are to describe and evaluate very high resolution DEMs of the studied rock glacier, to determine its surface velocity fields at various temporal scales, to evaluate the spatial characteristics of the rock glacier kinematics and to discuss the geomorphological significance of the results.

After presenting the study site (Section 3) and the methods used (Section 4), we describe the datasets (Section 5) and their analysis (Section 6). These last two sections include: (i) an evaluation of the data quality by estimating the vertical and horizontal errors; (ii) the determination of the surface displacement fields at various temporal scales; and (iii) the analysis of the surface displacement changes through time and space. Finally, the significance of those results in terms of possible causes of the spatio-temporal changes of rock glacier surface displacements is discussed (Section 7).

Study Site

Main Characteristics

The Laurichard active rock glacier is located in a transitional climate zone between the northern French Alps, affected by an oceanic regime, and the southern French Alps, mostly subjected to the Mediterranean influence. Using regionalization results from the SAFRANmodel provided by [START_REF] Durand | Reanalysis of 44 Yr of Climate in the French Alps (1958-2002): Methodology, Model Validation, Climatology, and Trends for Air Temperature and Precipitation[END_REF], the mean zero-isotherm altitude (ZIA) in the Laurichard area over the 1958-2002 period can be estimated to be around 2450 m asl; whereas the mean annual sum of precipitation is around 1000 mm at 1800 m asl, and the mean snow cover duration is about 220 days per year [START_REF] Durand | Reanalysis of 47 Years of Climate in the French Alps (1958-2005): Climatology and Trends for Snow Cover[END_REF].

The debris that constitutes the rock glacier is predominantly of granitic lithology, falling down from a densely-fractured 600 m-high rock wall. This gravitational activity feeds the rock glacier with relatively large boulders, which constitute an open-work debris accumulation that largely covers the slopes around.

The Laurichard rock glacier has a rather simple morphology, with an 800 m-long and 100-200 m-wide tongue, flowing down between 2650 m a.s.l. at the rooting zone (foot of the rock wall) and 2450 m a.s.l. at the snout. The activity of the rock glacier is evidenced by: (i) the presence of transversal ridges and furrows; (ii) a steep and unstable latero-frontal talus; and (iii) locally, by the looseness and the instability of the open-work active layer.

Previous Studies and Monitoring Activity

The Combe de Laurichard was chosen by Francou [START_REF] Francou | Chutes de pierres et éboulisation dans les parois de l'étage périglaciaire[END_REF] for achieving different geomorphological studies of scree slopes dynamics. Among other activities, Francou originally marked more than 40 blocks, among which 15 have been regularly measured since 1984 up to present (yearly since 1999, every 1-3 years before). Initially, blocks were positioned using a tacheometer [START_REF] Francou | 10 year surficial velocities on a rock glacier (Laurichard, French Alps)[END_REF], then a total station was used from 1992-2013, and since 2012, a differential dual-frequency GPS has been used in real-time kinematic (RTK) mode. This still ongoing long-term survey has allowed quantifying the annual surface velocity of the Laurichard rock glacier along its main flow line for more than three decades [START_REF] Bodin | Two decades of responses (1986-2006) to climate by the Laurichard rock glacier, French Alps[END_REF][START_REF] Bodin | Mountain permafrost and associated geomorphological processes: recent changes in the French Alps[END_REF] with an overall centimetric accuracy.

Comparison between the 2005 TLS datasets and a digital elevation model (DEM) representing the topography in 1975 (from digitizing the 10-m contours of the 1:25,000 IGN topographic map) shows that a large part of the rock glacier movement occurs on its right-hand side [START_REF] Bodin | High-Resolution DEM Extraction from Terrestrial LIDAR Topometry and Surface Kinematics of the Creeping Alpine Permafrost: The Laurichard Rock Glacier Case Study (Southern French Alps)[END_REF]. To account for this specific pattern of the rock glacier flow, new marked blocks were therefore added in 2008 to the original network of points (Figure 1).

Data and Methods

Datasets Collection

Since 2005, we conducted several measurement campaigns on the rock glacier aiming at generating DEMs that fulfill the following requirements: (i) a wide spatial coverage, at least the lower half of the rock glacier; (ii) a spatial resolution finer than the coarser blocks at the surface, the size of which is roughly estimated to be greater than 50 cm; and (iii) a quality that will allow at least an annual inter-comparison of the DEMs.

For that purpose, we performed several TLS and SFM-MVS campaigns. We also benefited from an ALS survey performed in 2012. After post-processing the initial datasets (see the following section), six DEMs have been generated at a resolution ranging from 5 cm-50 cm.

All DEMs were referenced in the RGF93 grid coordinate system based on the Lambert 93 conformal conic projection and the RAF09 geoid model for elevation reference (standard coordinate reference system of the French IGN, Institut National de l'Information Géographique et Forestière).

TLS Datasets

For the TLS surveys, we used an Optech R Ilris 3D scanner (Teledyne Optech,Vaughan, ON, Canada) (in 2005, 2006 and 2011) and an Optech R Ilris 3D-LR (long range) scanner (Teledyne Optech, Vaughan, ON, Canada) (in 2013). Details on the measurement capacities of the devices are described in Table 1.

For each of the four TLS surveys, between one and three stations were used to minimize the areas not viewed by the laser.

Between 5 and 12 artificial (35 cm-diameter spheres fixed on stable and moving rocks) and/or natural targets were measured during each campaign, using dual-frequency GNSS receivers, in differential mode (DGPS), allowing a relative accuracy of a few centimeters. The flight was performed by Sintégra R on 16 and 17 August 2012, with a Riegl LMS Q680i laser scanner on board a Cessna 206, flying at a nominal elevation above the ground of 610 m. The sky was clear with limited wind (<10 kts). The point density after filtering was 7 pts/m 2 . The embedded dual-frequency GNSS and inertial motion unit, coupled with a GNSS base station allow a positional accuracy of 5 cm and an angular accuracy of 0.01 • of the laser scanner position at each time during the flight.

Sintégra performed the post-processing (point cloud modeling, filtering, DEM generation and orientation), as well as the evaluation of the dataset by comparing to in situ GNSS measurements (45 points on a flat artificial surface, yielding a vertical accuracy (RMSE) of the 3D point cloud within 2 cm) [START_REF] Méchin | Acquisition D'une Couverture Topographique par MéThode Laser[END_REF].

SFM-MVS Datasets

On 16 September 2015, 27 pictures were taken from the Pyramide de Laurichard, i.e., facing the rock glacier (cf. Figure 1), with a Nikon D800 reflex camera, equipped with a 35-mm fixed-focal lens. The weather was cloudy, with mostly diffuse solar radiation, allowing a rather constant lighting of the terrain and homogeneous white balance distribution during the survey.

The initial raw (.nef Nikon file format) images were converted to .tiff format using XN Converter freeware. The images were processed in the commercial Photoscan Agisoft R software, with the following procedure: (1) selection of images using a quality index (7 images removed); (2) camera alignment and markers retrieval; (3) dense cloud processing; (4) export DEM at a 10-cm pixel resolution. Georeferencing of the 3D model is based on DGPS measurements of the artificial targets (spheres) and additional natural targets (blocks) visually detected on the 2012 ALS-derived DEM.

Datasets Processing

Co-Registration of the DEMs

Once each DEM had been individually georeferenced, they were all co-registered on stable areas with the 2005 DEM used as a reference. This was performed using a standard best-fit alignment procedure in InnovMetric R Polyworks software (InnovMetric Software Inc., Quebec, QC, Canada), which minimizes by a least-square method the distance between the two 3D models. The accuracy of the co-registration is assessed by the statistics of the error distribution (Table 2). In order to compute 3D displacements of the rock glacier between each pair of DEMs (2006-2005, 2011-2006, 2012-2011 and 2015-2012; note that the 2013 DEM appeared as too noisy and unfavorable for image correlation), the elevation grids were processed in the IMCORRmodule within SAGA-GIS [START_REF] Conrad | System for Automated Geoscientific Analyses (SAGA) v. 2.1.4[END_REF][START_REF] Scambos | Application of image cross-correlation to the measurement of glacier velocity using satellite image data[END_REF]. The algorithm retrieves pixels pattern between two georeferenced images and produces shapefiles (points and lines) containing the 2D and 3D displacements. Three parameters can be set: the search chip size, the pattern size and the grid sampling interval.

All the original DEMs were resampled at a 10-cm resolution and set in a common 2201 × 2951 pixels grid system. Several parameter combinations were tested, the best scores being achieved for the following parameter values: search chip: 64 or 128 pixels; pattern: 16 or 32 pixels; grid interval: 3 or 4.

The surface displacement datasets were manually cleaned due to border effects between grids with different local spatial coverage and artifacts (e.g., related to random local similarity of the coarse blocky surface) and then exported into a point-shapefile. This latter dataset was then interpolated (inverse distance weighting method, with a search radius of 10 m and a Gaussian weighting) into a 2-m resolution grid.

Comparison of Computed Surface Displacements with Ground-Based Measurements

As mentioned previously, 35 blocks located along three lines crossing the rock glacier (Figure 1) have been annually measured by the Parc national des Ecrins since 1999 [START_REF] Bodin | Two decades of responses (1986-2006) to climate by the Laurichard rock glacier, French Alps[END_REF]. Interannual displacements are therefore available for those points and can be compared, in the three dimensions, to the surface displacements (distance and direction) computed from images' correlation.

Results

Evaluation of the Vertical Accuracy

To assess the vertical accuracy of the datasets, we performed a pixel-to-pixel difference of the successive pairs of DEMs on stable areas (Table 3). For the 600-800 thousands points compared, the absolute maximal elevation difference is below 2.5 m, with a standard-deviation (SD) and a root-mean-squared error (RMS) within 15 and 27 cm. As expected, the distribution of the DEMs' difference is centered (the average equals zero), except for the 2015-2013 DEMs comparison, for which the 2015 model is locally suffering from broad-scale disturbances probably due to an uncorrected radial distortion inherent to an insufficient calibration in the SFM-MVS process [START_REF] James | Sequential digital elevation models of active lava flows from ground-based stereo time-lapse imagery[END_REF].

Evaluation of the Displacement Quantification

A preliminary visual inspection of the shaded-relief DTMs made from the 10-cm (2005,2006,2011,2013,2015) and 50-cm (2012) DEMs (cf. Supplementary Material) clearly shows the coherent patterns of the rock glacier deformation: stable areas outside the rock glacier, main central flow line along the maximal slope with reduced displacements towards the sides, higher velocity in the steepest upper part of the rock glacier and a bend of the movement towards its right orographic side.

Mapping the surface displacements derived from image correlation and from ground measurements shows a good agreement between the two approaches (Figure 2). The direction of the vectors is satisfactorily reproduced, as well as the their amplitude with an RMS error of 0.16 m/yr (Figure 3).

Interpretation of the Results

Spatio-Temporal Variation of the Rock Glacier Velocity

Though varying in amplitude, the spatial distribution of the movement remained rather homogeneous in space, with a very fast zone (>1 dm/yr) in the upper surveyed area and a very slow area (a few m/yr or less) on the left flank of the tongue (Figure 4). The orographic-right bend is also visible at every period, as well as marked border effects with a clear velocity lowering from the center flow line to the lateral sides. A localized lack of displacement vectors on the 2006-2011 datasets is due to a shallow dry debris-flow that took place over the steepest part of the rock glacier. This drastic morphological change, visible on the shaded-relief DTMs, impeded pattern retrieval by the image correlation algorithm. The mean ground-measured velocity (bold dashed line in Figure 5) and the velocity fields computed from image correlation over the whole rock glacier surface (red bold line in Figure 5) both display a close mean value, suggesting that the design of the annual topographic survey is well representative of the average rock glacier movement. Deviations of ±5-15% are also observable between both curves. Though this is within the error range of the velocity derived from image correlation, it can nevertheless be assumed that this could come from the heterogeneous velocity patterns: temporary acceleration or deceleration of specific areas of the rock glacier would be captured by the glacier-wide measurements derived from DEM correlation, but not by the network of blocks surveyed. The spread (standard deviation) of the image correlation-derived values also shows the overall large spatial variability of the velocity, especially changing between the rapid central flow line and the slower sides.

The longitudinal velocity profile (Figure 6

Deformation Patterns

The deformation pattern is here quantified by the relative spatial variation of velocity along the longitudinal profile (Figure 7), which mostly allows one to identify the compressive and extensive sectors along the profile. A compression area is visible near the front (horizontal distance > 200 m), due to the decrease in velocity just uphill of the front and the consequent buttressing effect. Higher up along the profile, a lateral shift of the compressive-extensive zone seems also observable between the two periods considered in our measurements, amounting to approximately 10 m, in coherence with the longitudinal translation of the ridges and furrows. Some significant changes, between the two periods, of the longitudinal velocity variation can be seen around 115 and 125 m (-20%) and around 130 and 140 m (+15%).

Discussion

Representativeness of the Estimated Surface Velocity

When statistically comparing the variability of the velocity measured from topographic surveys of marked blocks and the one of the velocity derived from image correlation of very high-resolution DEMs, we notice that: (i) the mean values of spatial variability from both methods are very close (0.21 m/yr for velocity derived from DEMs, 0.23 m/yr for velocity derived from topographic surveys), which means that the network of blocks is well representative of the whole behavior of the rock glacier; (ii) the temporal variability is also quite similar; therefore, even if the time interval between DEMs can be longer than the annual time span between two topographic surveys, a similar temporal signal is captured by both methods; (iii) the non-stationarity of the two methods is well identified: the mean velocity derived from DEMs increases, as well as its variability.

Topographic Control and Rheological Properties

No clear control of slope (derived from the 2012 DEM resampled to 2 m) on surface velocity can be detected when considering the entire profile (Figure 8). In order to account for the possible effect of surface roughness (many blocks are bigger than the 2-m DEM pixel size) and of the interannual variation of velocity, we also compared mean annual surface velocity measured with GPS between 2003 and 2016 with slopes derived from a 10-m DEM: similarly, no clear relationship appeared. Nevertheless, when considering only the 50-150-m section of the longitudinal profile (plain line section in Figure 7), a significant relationship is found for slope higher than 20 • (Figure 8). This specific pattern may be related to the specific rheological characteristics of the ice-debris mixture in a section of the profile that is mostly affected by creep, without acceleration and deceleration related to slope variation along the main flow line.

Limits and Perspectives of the Methods

Very high resolution DEMs, obtained from laser scanning or photogrammetry, have already been used for rock glacier studies. For instance, [START_REF] Buchli | Kinematic investigations on the Furggwanghorn Rock Glacier, Switzerland[END_REF] used multi-temporal LiDAR DEMs to derive velocity fields on a very fast rock glacier (surface velocities 2-7-times higher than the Laurichard case), whereas [START_REF] Kenner | Monitoring of high alpine mass movements combining laser scanning with digital airborne photogrammetry[END_REF] compared TLS-derived DEMs to airborne orthophoto on a variety of moving landforms, including rock glaciers. Nevertheless, our study is the first to assess the accuracy of LiDAR-and SFM-MVS-derived DEMs for quantifying multi-temporal surface displacements of rock glacier with 'normal' velocities (<2 m/yr).

We showed that interannual surface displacement fields of the Laurichard rock glacier can be derived in a rather straightforward way using high resolution DEMs. Though it is less precise than the geodetic survey of marked blocks, it is a complementary approach that allows describing precisely the spatial characteristics of the displacements. The combination of various sources of DEMs (terrestrial LiDAR, aerial LiDAR, SFM-MVS) is also possible, though it may complicate the selection of correct parameter values for the best retrieval of surface movement with the image correlation algorithm. One limit of this approach is the incomplete spatial coverage of terrestrial-based surveys, which leads to a lack of data in masked/hidden areas due to the limited view from the terrain. A good alternative could be, as proposed by [START_REF] Dall'asta | Unmanned Aerial Systems and DSM matching for rock glacier monitoring[END_REF], to use an unmanned aerial vehicle to produce a high resolution DEM that covers most of the study area.

Conclusions

By evaluating the very high resolution DEMs and surface displacement maps derived from those datasets, we can conclude that the various methods used to acquire the initial data (TSL, ALS, SFM-MVS) and the image correlation employed for further deriving displacement fields from DEMs are appropriate for investigating the multi-temporal kinematics of an active rock glacier. The accuracy we achieved for the DEMs, with horizontal and vertical errors around 0.15-0.3 m (excluding the ALS survey) and for the velocity fields (error of 0.16 m/yr), would also even allow one to quantify changes over a time interval shorter than a year, for instance to quantify displacement during summer. If both methods show similar spatial and temporal variability and are thus adapted for monitoring the pluri-annual velocity changes, temporary and localized changes in velocity, not captured by the topographic surveys of the blocks, may also occur as suggested by the differences between mean values derived from both methods. Therefore, in order to investigate the seasonal behavior of rock glaciers and to detect localized changes, we suggest to perform, every 1-3 years, TLS, ALS or terrestrial/drone SFM-MVS surveys in order to produce 10-50-cm DEMs, in addition to the annual topographic survey.

Figure 1 .

 1 Figure 1. Map of the study area with the location of: (1) the targets (spheres with a diameter of 30 cm, anchored on stable terrain) used for georeferencing; (2) the annually surveyed marked blocks (DGPS); and (3) the cameras for the 2015 survey (symbols are oriented with the optical axis of the lens). The lower-left inset shows the general location of the study area in the French Alps (source: Open Street Map).

Figure 2 .

 2 Figure 2. Surface velocity (m/yr) on the Laurichard rock glacier for the period 2011-2012 quantified by image correlation on DEMs (red lines) and by ground topographic (total station and DGPS) measurements of marked blocks (black arrows).

Figure 3 .

 3 Figure 3. Surface velocity derived from image correlation (computed from the values of the 10 closest points around the actual position of the marked blocks) versus 3D surface velocity measured on marked blocks.

Figure 4 .

 4 Figure 4. Maps of the surface velocity of the Laurichard rock glacier at three different time periods, as derived from image correlation using HR DEMs.

  ) shows a generalized progressive increase in the surface displacements, between 2005-2006 and 2012-2015. The mean velocity along the profile changed from 0.82-1.80 m/yr, which corresponds to a 45% increase over the entire period (respectively 6, 27 and 21% for 2005-2006/2006-2011, 2006-2011/2011-2012, 2011-2012/2012-2015). This acceleration is not homogeneously distributed along the profile and varies through time: it remained almost constant (around 40-50%) along the steepest part (horizontal distance: from 50-150 m, with partial data gaps of 25 and 80 m). On the lowest part (horizontal distance: above 150 m), the velocity remained rather constant (interannual velocity variation roughly ranging from 30--30%) over the first two periods (2005-2006/2006-2011, 2006-2011/2011-2012), whereas it clearly increased between 2011-2012 and 2012-2015.

Figure 5 .

 5 Figure 5. Evolution of the Laurichard rock glacier mean velocity since the early 2000s, derived from ground topographic measurements (total station and DGPS; bold black dashed line) and from image (DEM) correlation (red lines).

Figure 6 .

 6 Figure 6. Variation of velocity along the longitudinal axis of the Laurichard rock glacier between 2005 and 2015, derived from image correlation. The upper panel shows the velocity profile for each period, and the lower panel depicts the variations of velocity between the different periods.

Figure 7 .

 7 Figure 7. Relative spatial variation of velocity along the longitudinal axis of the Laurichard rock glacier for the slowest period (2005-2006, in black) and the fastest period (2012-2015, in red). Positive values indicate an increasing velocity towards downslope (extensive deformation), whereas negative values indicate a slowing down (compressive deformation).The plain line is used in the central section of the profile, which is referred to later on (see Figure8).

Figure 8 .

 8 Figure 8. Velocity against slope along the longitudinal profile (circle), with values for the central section (horizontal distance between 50 and 150 m; black dot). The red line is the regression fitted to the values for the central section for a slope angle above 20 • .

Table 1 .

 1 Technical overview of the terrestrial laser scanners used in the study (from the Optech R factsheets).

		Ilris 3D	Ilris 3D-LR
	Range (80% reflectivity)	1200 m	3000 m
	Range (10% reflectivity)	400 m	1300 m
	Accuracy on distance	7 mm @100 m	7 mm @ 100 m
	Angular accuracy	80 µrad	80 µrad
	Position accuracy	8 mm @ 100 m	8 mm @ 100 m
	Laser footprint	29 mm @ 100 m 27 mm @ 100 m
	Laser beam divergence	0.17 mrad	0.25 mrad
	Wavelength	1500 nm	1064 nm
	Acquisition frequency	2 kHz	10 kHz
	Field of view	40 • × 40 •	40 • × 40 •
	4.1.2. ALS Datasets

Table 2 .

 2 Accuracy of the datasets, with absolute values based on the manufacturer factsheet (TLS) or software processing results (SFM-MVS with Agisoft Photoscan), and relative values based on the co-registration error of each DEM on 2005 DEM. * Data according to the manufacturer factsheets.

		2005	2006	2011	2012	2013	2015
	Type of survey	TLS	TLS	TLS	ALS	TLS	SFM-MVS
	Date of survey	18 September	25 August	13 September	16/17 August	10 October	16 September
	Sensor	Ilris-3D	Ilris-3D	Ilris-3D	Riegl LMS Q680i Ilris-3D LR	Nikon D800
	Position accuracy *	8 mm @ 100 m 8 mm @ 100 m 8 mm @ 100 m 20 mm @ 250 m	7 mm @ 100 m 68 mm (6 GCPs)
	Co-registration
	error on 2005 point	-	0.102	0.030	0.003	-0.053	0.144
	cloud: mean (m)
	Co-registration
	error on 2005 point	-	0.301	0.289	0.513	0.213	0.148
	cloud: SD (m)
	Number of non-NApoints	-	332,831	287,860	14,029	100,992	903,269
	compared
	4.2.2. Image Correlation for Surface Displacement Measurement

Table 3 .

 3 Vertical differences (in meters) between successive pairs of DEMs in stable areas around the rock glacier (sub-vertical rock face, flat Quaternary moraine and gently inclined fossil rock glacier).

		2006-2005 2011-2006 2012-2011 2013-2012 2015-2013
	Min	-2.43	-2.4	-2.47	-2.02	-1.82
	Max	2.1	2.29	2.48	1.4	2.27
	Median	0.03	-0.03	-0.00	0.03	0.19
	Mean	0.03	-0.01	-0.04	0.04	0.18
	SD	0.25	0.21	0.23	0.14	0.20
	RMSE	0.25	0.21	0.24	0.15	0.27
	non-NA	732,000	766,023	832,995	642,490	635,736

Remote Sens. 2018, 10, 547; doi:10.3390/rs10040547 www.mdpi.com/journal/remotesensing

Acknowledgments:

We acknowledge the Parc national des Ecrins for operating the topographic survey of the rock glacier since 1999. We are also grateful for the many people who helped in the field and for the SAJF (UMS 2925 CNRS/UGA), who provided us with the 50-cm ALS DEM and sometimes hosts us. The Ilris3D-LR laser scanner used in the study has been funded by the LabEx OSUG@2020 (Investissements d'avenir-ANR10 LABX56).

Supplementary Materials:

The following are available online at www.mdpi.com/2072-4292/10/4/547/s1. The six 0.1-m pixel-size DEMs used in this study are available as Supplementary Material.

Author Contributions: X.B. coordinated the study, the data processing and interpretation and the writing of the manuscript. E.T. performed the comparison analysis of the spatial-temporal variability of both methods. O.S. processed the point clouds and performed the co-registration. A.R. helped in the TLS data acquisition, the interpretation of the results and edited the manuscript. S.J. helped with the acquisition and the processing of the TLS datasets.

Conflicts of Interest:

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results.


Adrien Michez
email: adrien.michez@uliege.be

Philippe Lejeune
email: p.lejeune@uliege.be

Sébastien Bauwens
email: bauwens.sebastien@gmail.com

Andriamandroso Andriamasinoro
email: herinaina.andriamandroso@yncrea.fr

Lalaina Herinaina

Yannick Blaise
email: blaiseyannick@hotmail.com

Eloy Castro Muñoz

Frédéric Lebeau
email: f.lebeau@ulg.ac.be

Jérôme Bindelle
email: jerome.bindelle@uliege.be

remote sensing

Keywords: unmanned aerial vehicles, unmanned aerial systems, drone, precision grazing, pasture biomass modeling, sward height, pasture height

The tools available to farmers to manage grazed pastures and adjust forage demand to grass growth are generally rather static. Unmanned aerial systems (UASs) are interesting versatile tools that can provide relevant 3D information, such as sward height (3D structure), or even describe the physical condition of pastures through the use of spectral information. This study aimed to evaluate the potential of UAS to characterize a pasture's sward height and above-ground biomass at a very fine spatial scale. The pasture height provided by UAS products showed good agreement (R 2 = 0.62) with a reference terrestrial light detection and ranging (LiDAR) dataset. We tested the ability of UAS imagery to model pasture biomass based on three different combinations: UAS sward height, UAS sward multispectral reflectance/vegetation indices, and a combination of both UAS data types. The mixed approach combining the UAS sward height and spectral data performed the best (adj. R 2 = 0.49). This approach reached a quality comparable to that of more conventional non-destructive on-field pasture biomass monitoring tools. As all of the UAS variables used in the model fitting process were extracted from spatial information (raster data), a high spatial resolution map of pasture biomass was derived based on the best fitted model. A sward height differences map was also derived from UAS-based sward height maps before and after grazing. Our results demonstrate the potential of UAS imagery as a tool for precision grazing study applications. The UAS approach to height and biomass monitoring was revealed to be a potential alternative to the widely used but time-consuming field approaches. While reaching a similar level of accuracy to the conventional field sampling approach, the UAS approach provides wall-to-wall pasture characterization through very high spatial resolution maps, opening up a new area of research for precision grazing.

Introduction

Grasslands cover 26% of the world's total land area and 69% of the agricultural area and are the least expensive way to feed ruminants, such as cattle reared for milk or meat production [START_REF] O'mara | The role of grasslands in food security and climate change[END_REF]. In the European Union, permanent pastures and meadows occupy 34% of the utilized agricultural area. A significant share of the cultivated fields include temporary pastures in their rotation [START_REF]Eurostat Farm Structure Statistics-Statistics Explained[END_REF]. Grazed pastures benefit the sustainability of dairy production in multiple ways, such as lower feeding costs [START_REF] Dillon | Optimising financial return from grazing in temperate pastures[END_REF], higher animal welfare and a lower occurrence of lameness and mastitis, good public image, and increased milk quality [START_REF] Burow | Effect of grazing on the cow welfare of dairy herds evaluated by a multidimensional welfare index[END_REF]. In addition, grasslands have an important role in the provision of social and environmental services [START_REF] Boval | The importance of grasslands for animal production and other functions: A review on management and methodological progress in the tropics[END_REF]. The selection of an adequate choice of diverse grass and legume species and varieties with adequate management could also support the growth of a wider range of micro-fauna and crop auxiliaries. Finally, pastures could play a significant role in trapping atmospheric CO 2 through soil carbon sequestration [START_REF] Dumortier | Methane balance of an intensively grazed pasture and estimation of the enteric methane emissions from cattle[END_REF].

Managing pastures through direct grazing is not an easy task, as it requires combining actions to allow the dynamic growth of multispecies grassland vegetation while managing animals that consume it selectively and non-homogeneously. The plant and the animal components of the grazed pasture system have mutual influences in addition to the influences of environmental factors. Both components also have different biological dynamics. Unfortunately, the tools available to farmers to manage grazed pastures and adjust forage demand to grass growth are generally rather static, for example, the use of stocking rates, fences, grazing pressure, resting times, locations of water troughs, locations of salt licks, fertilization, and irrigation [START_REF] Holechek | Principles and Practices[END_REF].

The application of precision livestock farming to grazing has the potential to promote more dynamic management of ruminants on grasslands by switching from management at the herd level to the management of each individual on the pasture. In the past decade, several works have demonstrated how the development in sensors and information technology can be used to improve the monitoring of grazing animals, especially cattle [START_REF] Andriamandroso | A review on the use of sensors to monitor cattle jaw movements and behavior when grazing[END_REF][START_REF] Debauche | Web-based cattle behavior service for researchers based on the smartphone inertial central[END_REF][START_REF] Laca | Precision livestock production: Tools and concepts[END_REF][START_REF] Larson-Praplan | Spatial and temporal domains of scale of grazing cattle[END_REF]. Monitoring the grazing process to the finest scale, namely at the consumption level, in space and time is now possible [START_REF] Andriamandroso | A review on the use of sensors to monitor cattle jaw movements and behavior when grazing[END_REF], and the development of such technology as a research tool is likely to enter farms progressively as a complement to virtual fences [START_REF] French | Development and adoption of new technologies to increase the efficiency and sustainability of pasture-based systems[END_REF]. In addition to the monitoring of individual grazing animals, improvements in the monitoring of the grazed vegetation are also called for. French et al. [START_REF] French | Development and adoption of new technologies to increase the efficiency and sustainability of pasture-based systems[END_REF] presented a rising plate meter that automatically communicates the pasture height and biomass data to a geographic information system (GIS). Other techniques, such as satellite-based remote-sensing [START_REF] Schellberg | Frontiers and perspectives on research strategies in grassland technology[END_REF] and on-field 2-D and 3D-cameras, have been used to measure photosynthetic activity and forage height [START_REF] Andriamandroso | Differentiating pre-and post-grazing pasture heights using a 3D camera: A prospective approach[END_REF] or to identify weeds in grasslands [START_REF] Van Evert | Real-time vision-based detection of Rumex obtusifolius in grassland[END_REF]. Handcock et al. [START_REF] Handcock | Monitoring animal behaviour and environmental interactions using wireless sensor networks, GPS collars and satellite remote sensing[END_REF] proposed combining satellite remote sensing and collar-based georeferenced information to characterize animal behavior and environmental interactions. Even though the accessibility and the spatial resolution of satellite imagery are continuously improving, generally, they are too coarse to describe pastures at the scale of individual animals on the pasture. On the other hand, field-based remote sensing methods, such as terrestrial light detection and ranging (LiDAR), allow for the fine 3D characterization of vegetation structure [START_REF] Tilly | Multitemporal crop surface models: Accurate plant height measurement and biomass estimation with terrestrial laser scanning in paddy rice[END_REF]. LiDAR technology was also successfully combined with a vegetation index (normalized difference vegetation index, NDVI) by Schaeffer and Lamb [START_REF] Schaefer | A combination of plant NDVI and LiDAR measurements improve the estimation of pasture biomass in tall fescue (Festuca arundinacea var. Fletcher)[END_REF] to estimate the standing biomass on tall fescue pastures. Nevertheless, terrestrial remote sensing remains spatially limited. At the interface of satellite-based remote sensing and on-field approaches, unmanned aerial systems (UASs) can cover substantial areas (>10 ha in a single flight), providing imagery at a very high spatial resolution (<0.1 m). Moreover, UASs are very versatile tools that can be deployed by the end user on-demand. Over the past decade, UASs have been widely used in various environmental applications including forest characterization [START_REF] Michez | Classification of riparian forest species and health condition using multi-temporal and hyperspatial imagery from unmanned aerial system[END_REF], wildlife census [START_REF] Vermeulen | Unmanned Aerial Survey of Elephants[END_REF], invasive species mapping [START_REF] Michez | Mapping of riparian invasive species with supervised classification of Unmanned Aerial System (UAS) imagery[END_REF], and precision agriculture [START_REF] Zhang | The application of small unmanned aerial systems for precision agriculture: A review[END_REF]. In precision agriculture, UASs have been used for quantitative crop monitoring [START_REF] Michez | How Far Can Consumer-Grade UAV RGB Imagery Describe Crop Production? A 3D and Multitemporal Modeling Approach Applied to Zea mays[END_REF], weed detection [START_REF] Lopez-Granados | Weed detection for site-specific weed management: Mapping and real-time approaches[END_REF], and plant phenotyping [START_REF] Yang | Unmanned Aerial Vehicle Remote Sensing for Field-Based Crop Phenotyping: Current Status and Perspectives[END_REF]. They have also been used for grazing applications in research, including biomass monitoring [START_REF] Lee | Mapping Herbage Biomass on a Hill Pasture using a Digital Camera with an Unmanned Aerial Vehicle System[END_REF] and adventive species detection [START_REF] Lee | Estimating the spatial distribution of Rumex acetosella L. on hill pasture using UAV monitoring system and digital camera[END_REF]. The success of UASs is linked to the development of structure from motion (SFM) photogrammetry [START_REF] Westoby | Structure-from-Motion" photogrammetry: A low-cost, effective tool for geoscience applications[END_REF], which allows 3D data (digital surface model in case of aerial survey) to be derived from overlapping images acquired with consumer-grade digital cameras. In the agricultural field, SFM photogrammetry can be used to derive the crop surface/height model (CSM/CHM) from UAS imagery at very low cost compared to LiDAR. In the context of precision grazing, pioneer studies [START_REF] Bareth | Feasibility study of using non-calibrated UAV-based RGB imagery for grassland monitoring: Case study at the Rengen Long-term Grassland Experiment (RGE), Germany[END_REF][START_REF] Viljanen | A Novel Machine Learning Method for Estimating Biomass of Grass Swards Using a Photogrammetric Canopy Height Model, Images and Vegetation Indices Captured by a Drone[END_REF] have demonstrated the quality of height and biomass information derived from the UAS. However, these previous studies have mostly focused on ungrazed, experimental pastures while leaving open the question of monitoring grazing through the use of a multi-temporal approach. Such approaches are needed to develop innovative tools to provide quick responses to allow the adaptation of grazing schedules or other decision-making tools for farmers.

Pastures should be consumed by animals when they display a specific structure that is principally defined by height and density, ensuring the removal of only grass leaves during bites [START_REF] Da Trindade | Forage allowance as a target of grazing management: Implications on grazing time and forage searching[END_REF]. Pasture height information is very important for farmers to maximize the intake rate of animals. However, such measurements are still taken using low frequency and time-consuming instruments, notably sward-sticks or rising plate meters [START_REF] Andriamandroso | Differentiating pre-and post-grazing pasture heights using a 3D camera: A prospective approach[END_REF]. The use of UAS imagery could be a solution to provide higher frequency, higher spatial resolution, and more precise information about the biomass productivity of a pasture.

In this context, this study aimed to evaluate the potential of UAS as a tool for the characterization of pasture 3D structure (sward height) and above-ground biomass at a very fine spatial scale. More specifically, it evaluated the suitability of UAS for monitoring sward height differences and pasture biomass. Based on the hypothesis that UAS products (e.g., sward height model and orthophotos mosaic) can reflect the pasture biomass and sward height multi-temporal differences, this approach investigated the potential of UASs to monitor grazing at the feeding station level.

Materials and Methods

Study Site and Field Sampling

The study site consisted of experimental Lolium perenne (ryegrass) and Trifolium repens (clover) based pastures (Figure 1) located in Gembloux (Belgium). Each individual plot represented an area of ca. 4000 m 2 . The study took place from 9 May to 12 May 2017. All pastures had been deferred from grazing since October 2016 (seven months), and grass was allowed to regrow ungrazed until the beginning of the experiment. Six dry red-pied cows were introduced to the first pasture (see Figure 1, plot P1) for two consecutive days (9 May to 10 May) and were moved to the next test pasture (P2) for the last two days (11 May to 12 May). Forty above-ground biomass samples were taken for dry-weighing from both pastures (P1 and P2) within a square of 0.09 m 2 centered on the marked point (four squares and center) along the sampling grid (Figure 1) before cows were set to graze. No distinction between the different species was made in the data collection. The center of the biomass samples was precisely geo-referenced using an RTK GPS (±0.04 m mean XYZ accuracy).

UAS Imagery Acquisition and Processing

An octocopter drone (X frame type) equipped with an off-the-shelf high resolution (20 Mpx) RGB camera (Sony RX 100, 8.8 mm focal length) and a multispectral camera (Parrot Sequoia) was used. The multispectral camera had a lower spatial resolution (1.2 Mpx), but provided fine multispectral information, covering green (550 nm), red (660 nm), near infrared (735 nm), and red-edge (790 nm) wavelengths. The full width at half maximum was 40 nm for red, green, and near infrared, and 10 nm for the red-edge wavelength. The flight height was set to 50 m above ground level, with a cruise speed of 5 m.s -1 and front and side overlaps above 80%. The UAS flight surveys were performed before (May 8) and after the cows grazed the pastures (15 May).

1, plot P1) for two consecutive days (9 May to 10 May) and were moved to the next test pasture (P2) for the last two days (11 May to 12 May). Forty above-ground biomass samples were taken for dry-weighing from both pastures (P1 and P2) within a square of 0.09 m² centered on the marked point (four squares and center) along the sampling grid (Figure 1) before cows were set to graze. No distinction between the different species was made in the data collection. The center of the biomass samples was precisely geo-referenced using an RTK GPS (±0.04 m mean XYZ accuracy). Eight ground control points (GCPs) consisting of 0.4 m white square plastic plates were placed in the surveyed area to ensure geometric calibration. The GCPs were georeferenced with an RTK GPS (±0.04 m mean XYZ accuracy).

The images provided by the high spatial resolution RGB sensor (Sony RX100 III) were processed with Agisoft Photoscan 1.3 to produce a high spatial resolution digital surface model (DSM, 0.025 m GSD). Agisoft Photoscan is widely used by the UAS scientific community to perform high resolution and reliable 3D models from UAS surveys (e.g., [START_REF] James | Optimising UAV topographic surveys processed with structure-from-motion: Ground control quality, quantity and bundle adjustment[END_REF][START_REF] Javernick | Modeling the topography of shallow braided rivers using Structure-from-Motion photogrammetry[END_REF][START_REF] Gonçalves | UAV photogrammetry for topographic monitoring of coastal areas[END_REF]). The sparse and dense clouds were processed using the "high" level of accuracy for parameters. The depth filtering strategy was set to "aggressive". The optimization process was based on the GCP location and was applied to all of the available parameters (except for the rolling shutter parameter). As Agisoft Photoscan 1.3 does not allow radiometric calibration workflow, we used Pix4D Suite 3.1 (Lausanne, Switzerland) to perform a photogrammetric 3D reconstruction of the multispectral imagery using the "Ag Multispectral" predefined workflow. Radiometric calibration of the multispectral camera was undertaken before each flight with a calibration target (Airinov). The images acquired by the multispectral sensor were used to derive reflectance maps (0.05 m ground sampling distance, GSD) of the four wavelengths (green, red, near-infrared, and red-edge). The four reflectance layers were combined to produce four straight-forward vegetation indices (0.05 m GSD) covering the different layers of the multispectral camera (Table 1). We did not include the spectral information provided by the visible camera in the analysis because of its lower spectral resolution.

Using the same approach as Michez et al. [START_REF] Michez | How Far Can Consumer-Grade UAV RGB Imagery Describe Crop Production? A 3D and Multitemporal Modeling Approach Applied to Zea mays[END_REF] with Zea mays field crops, a height model was derived by subtracting a LiDAR digital terrain model (DTM) from the high spatial resolution photogrammetric DSM in order to provide a raster of the sward height at the study site (0.025 m GSD). This height model will be further referred to as the sward height model (SHM). The LiDAR DTM used in this study was derived from a relatively low-density LiDAR survey (<1 point/m 2 ) acquired during the years 2013 and 2014 (see http://geoportail.wallonie.be/forfurtherinformation). As the pastures were installed before the LiDAR survey and have not been tilled since, this LiDAR DTM was considered to be relevant for describing the topography of the study site.

Remote Sens. 2019, 11, 473 5 of 14 Table 1. Vegetation indices computed from the reflectance layers provided by the multispectral camera.

Vegetation Index Formula Parameter Reference

Normalized Difference Vegetation Index (NDVI) (NIR -RED)/(NIR + RED) Photosynthetic activity, plant stress [START_REF] Rouse | Monitoring Vegetation Systems in the Great Plains with ERTS[END_REF] Normalized Difference Red Edge (NDRE) (NIR -REDEDGE)/(NIR + REDEDGE) Chlorophyll and nitrogen content [START_REF] Barnes | Coincident detection of crop water stress, nitrogen status and canopy density using ground based multispectral data[END_REF] Green NDVI (GNDVI) (NIR -GREEN)/(NIR + GREEN)

More sensitive to chlorophyll-a concentration, monitoring of plant stress [START_REF] Moges | Evaluation of green, red, and near infrared bands for predicting winter wheat biomass, nitrogen uptake, and final grain yield[END_REF] Green Ratio Vegetation Index (GRVI) NIR/GREEN Photosynthetic activity [START_REF] Sripada | Aerial Color Infrared Photography for Determining Early In-Season Nitrogen Requirements in Corn[END_REF] At the center of the field measurement locations, we extracted the UAS information (height and reflectance values) as the median values of the considered layers' pixels in a 0.15 m radius circular buffer. The UAS variables used for the analysis are listed in Table 2.

Validation of UAS Sward Height

Five static LiDAR laser scans were performed using a FARO Focus 3D 120 (FARO, 250 Technology Park Lake Mary, FL 32746, United States) to produce reference sward height data. The LiDAR scanner uses phase-shift-based LIDAR technology to measure the XYZ locations of objects surrounding the scanner. The FARO scanner has a maximum range of 120 m with an accuracy of 2 mm at 10 m. The Focus 3D 120 has a 360 • horizontal field of view and a 305 • vertical field of view. Due to the 55 • missing vertical field of view, a disc of 1.5 m of radius was not scanned under the scanner (the position height of the scanner was 1.3 m). The scanner was located in the center of five subsamples of the plots. The resolution of the scans was set to one-quarter of the full resolution. The different locations of the scanner were recorded with an RTK GPS (±0.04 m mean XYZ accuracy). The LiDAR scans were performed on May 7 before cows were set to graze.

The five generated point clouds were georeferenced using the FARO SCENE software (version 6.2.4). The LASTools suite (version 171124, http://rapidlasso.com/LAStools) was used to produce five digital surface models corresponding to a triangulation (0.01 m grid) of the highest points. In order to remove the noise, triangulation was performed over the 95th percentile points (Z values) within a XY grid of 0.01 m. The five DSM were subtracted to the above-referenced aerial LiDAR DTM to produce a LiDAR sward height model (SHM) at 0.01 m spatial resolution. This last step was performed to counter the low ground point density in the terrestrial laser scans due to occlusion by the pasture canopy.

The five grass LiDAR-derived SHM were compared to the UAS SHM (before grazing) using a 0.5 m grid within a circular area of interest (2 m < radius < 5 m) centered on the location of the scans (Figure 2). For every grid cell, a median value of the compared SHM was computed within a circular buffer of 0.1 m. A linear regression model was fitted to evaluate the correspondence between the sward height estimated with the UAS imagery and a reference SHM extracted with highly accurate LiDAR data:

Model 1: Ref. LiDAR sward height = f (UAS sward height).

(1)

Model 1 was a linear regression model fitted between the pasture height in the study area obtained from the LiDAR SHM and the UAS SHM (575 observations before cattle grazing). Considering the high quality of the LiDAR 3D products, the LiDAR pasture height was considered as a reference pasture height. The goodness of fit of Model 1 was considered as a proxy of the ability of the UAS height to describe the pasture height in further analyses, notably to produce sward height maps before and after grazing as well as sward height difference maps from multitemporal sward height models. A cross validation of Model 1 was performed using the repeated k-fold approach (k = 5; repetitions = 100).

Modeling Biomass of Pasture with UAS Imagery

The modeling of pasture biomass (before grazing) was investigated comparatively with different UAS products (3D and/or spectral):

Models 2, 3 and 4 evaluated the strictly 3D, spectral, and mixed (3D + spectral) UAS biomass modeling approaches, respectively (see Table 2 for a detailed list of variables). The UAS heights and spectral data averaged (median) within a circle of 0.15 m radius were compared to the pasture biomass values measured on the 0.09 m² georeferenced spots. We checked for non-linearity between the predictors and the field measured pasture biomass. To reduce the number of predictors in the fitted multivariate models (Models 3 and 4), stepwise selection (both directions) was performed using the Akaike information criterion [START_REF] Akaike | Maximum likelihood identification of Gaussian autoregressive moving average models[END_REF] implemented in the StepAIC tool from the "MASS" R package [START_REF] Venables | Modern Applied Statistics with S[END_REF]. Within the selected set of variables, only significant variables (as determined by p-value) were kept in the final model. To determine the pasture biomass model that presented the best performance, the relative importance of the selected variables was computed by a decomposition of the models' explained variance into non-negative contributions, following the approach of Lindeman et al. [START_REF] Lindeman | Introduction to Bivariate and Multivariate Analysis[END_REF]. This was implemented in the "relaimpo" package in R. For all of the UAS pasture biomass models, cross validation using the repeated k-fold approach was performed (k = 5; repetitions = 100).

Mapping Sward Height Differences and Biomass with UAS

As all of the UAS variables used in the model fitting process were extracted from spatial information (raster data), a high spatial resolution map of pasture height and pasture biomass was derived based on the previously fitted models. The multitemporal height maps were combined to

Models 2, 3 and 4 evaluated the strictly 3D, spectral, and mixed (3D + spectral) UAS biomass modeling approaches, respectively (see Table 2 for a detailed list of variables). The UAS heights and spectral data averaged (median) within a circle of 0.15 m radius were compared to the pasture biomass values measured on the 0.09 m 2 georeferenced spots. We checked for non-linearity between the predictors and the field measured pasture biomass. To reduce the number of predictors in the fitted multivariate models (Models 3 and 4), stepwise selection (both directions) was performed using the Akaike information criterion [START_REF] Akaike | Maximum likelihood identification of Gaussian autoregressive moving average models[END_REF] implemented in the StepAIC tool from the "MASS" R package [START_REF] Venables | Modern Applied Statistics with S[END_REF]. Within the selected set of variables, only significant variables (as determined by p-value) were kept in the final model. To determine the pasture biomass model that presented the best performance, the relative importance of the selected variables was computed by a decomposition of the models' explained variance into non-negative contributions, following the approach of Lindeman et al. [START_REF] Lindeman | Introduction to Bivariate and Multivariate Analysis[END_REF]. This was implemented in the "relaimpo" package in R. For all of the UAS pasture biomass models, cross validation using the repeated k-fold approach was performed (k = 5; repetitions = 100).

Mapping Sward Height Differences and Biomass with UAS

As all of the UAS variables used in the model fitting process were extracted from spatial information (raster data), a high spatial resolution map of pasture height and pasture biomass was derived based on the previously fitted models. The multitemporal height maps were combined to produce a high resolution map of sward height differences (height map 'before' MINUS height map 'after'). As field biomass sampling only occurred 'before grazing', the pasture biomass map was subsequently computed for this time window.

Results

Validation of UAS Sward Height

The sward height provided by UAS products agreed well (R 2 = 0.62, RMSE = 0.04 m) with the reference (LiDAR) sward height information (Figure 3). The results of the cross validation process showed a similar performance and can be found in Appendix A.

Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 14

Results

Validation of UAS Sward Height

The sward height provided by UAS products agreed well (R² = 0.62, RMSE = 0.04 m) with the reference (LiDAR) sward height information (Figure 3). The results of the cross validation process showed a similar performance and can be found in Appendix A.

Regarding our objectives, the accuracy of the sward height estimated with UAS imagery was consequently considered to be sufficient to produce sward height maps for the study site (before and after grazing) to derive a sward height differences map.

Modeling Biomass of Pasture with UAS Imagery

The quality of the fitted pasture biomass models (Figure 4) presented adjusted R² values ranging from 0.23 (Model 2) to 0.49 (Model 4) with a RMSE of ca. 0.1 kg/m². These results also show that the spectral information provided by the UAS imagery was meaningful, as the spectral biomass model (Model 3) displayed better performance than the UAS sward height biomass model (Model 2). The complementarity of spectral information and sward height information provided by UAS in the mixed model (Model 4) was also highlighted, as Model 4 had the best performance (adj. R² = 0.49, R² = 0.62). Despite the complementarity of the UAS height and the UAS spectral information, 45% of the variance of the mixed biomass model (Model 4) is associated with the UAS height (Table 3). The results of the cross validation (found in Appendix A) highlighted a similar performance for all of the pasture biomass UAS models. Regarding our objectives, the accuracy of the sward height estimated with UAS imagery was consequently considered to be sufficient to produce sward height maps for the study site (before and after grazing) to derive a sward height differences map.

Modeling Biomass of Pasture with UAS Imagery

The quality of the fitted pasture biomass models (Figure 4) presented adjusted R 2 values ranging from 0.23 (Model 2) to 0.49 (Model 4) with a RMSE of ca. 0.1 kg/m 2 . These results also show that the spectral information provided by the UAS imagery was meaningful, as the spectral biomass model (Model 3) displayed better performance than the UAS sward height biomass model (Model 2). The complementarity of spectral information and sward height information provided by UAS in the mixed model (Model 4) was also highlighted, as Model 4 had the best performance (adj. R 2 = 0.49, R 2 = 0.62). Despite the complementarity of the UAS height and the UAS spectral information, 45% of the variance of the mixed biomass model (Model 4) is associated with the UAS height (Table 3). The results of the cross validation (found in Appendix A) highlighted a similar performance for all of the pasture biomass UAS models.

Mapping Sward Height Differences and Biomass with UAS

The maps of the UAS sward height before and after grazing (Figure 5I,II) presented interesting spatial patterns with a clearly distinguishable area with higher values in the center of the study site in the 'before grazing' height map (Figure 5I). The combination of sward height maps 'before' and 'after' grazing allowed for the computation of a sward height differences map (Figure 5III). Considering the distribution represented in the density curve (Figure 5IV), it can be seen that most of the study site was grazed or trampled as the peak was associated with positive values. The area with a higher sward height valued in the 'before' map (Figure 5I), located in the center of the area, was more intensively grazed or trampled.

Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 14

The maps of the UAS sward height before and after grazing (Figures 5I,II) presented interesting spatial patterns with a clearly distinguishable area with higher values in the center of the study site in the 'before grazing' height map (Figure 5I). The combination of sward height maps 'before' and 'after' grazing allowed for the computation of a sward height differences map (Figure 5III). Considering the distribution represented in the density curve (Figure 5IV), it can be seen that most of the study site was grazed or trampled as the peak was associated with positive values. The area with a higher sward height valued in the 'before' map (Figure 5I), located in the center of the area, was more intensively grazed or trampled. Figure 6 allowed us to visualize the spatial pattern of the pasture biomass and its heterogeneity at a very fine scale. This map was produced through the use of Model 4 applied to the associated UAS variables (i.e., UAS sward height, GRVI, GNDVI, and NDRE). The northwestern part of the study site presented lower biomass values. Some abnormal values were shown on the empirical cumulative distribution function (ECDF) plot (Figure 6II) that were well above the usual averages observed for such pastures or were even negative. Nevertheless, these values can be considered as noise when considering that the 5th and 95th percentiles of the predicted (and mapped) pasture biomass values were 0.04 kg/m² and 0.36 kg/m², respectively. 6II) that were well above the usual averages observed for such pastures or were even negative. Nevertheless, these values can be considered as noise when considering that the 5th and 95th percentiles of the predicted (and mapped) pasture biomass values were 0.04 kg/m 2 and 0.36 kg/m 2 , respectively.

Discussion

Modeling Biomass of Pasture with UAS Imagery

The performance (adj. R² = 0.49) of the best performing UAS-based biomass model (Model 4) was similar to the state-of-art field pasture biomass non-destructive monitoring tools for swards-turf grasses such as the rising plate meter [START_REF] Harmoney | Determination of pasture biomass using four indirect methods[END_REF][START_REF] Laca | Comparison of weight estimate and rising-plate meter methods to measure herbage mass of a mountain meadow[END_REF]. There are still very few references for the application of UAS imagery products to pasture biomass modeling on such intensively managed pastures with rather low pasture heights. Lee et al. [START_REF] Lee | Mapping Herbage Biomass on a Hill Pasture using a Digital Camera with an Unmanned Aerial Vehicle System[END_REF] achieved a better model performance (R² = 0.77) than this study, but their work was performed on a multisite dataset in Korea, which had a higher diversity and more contrasting biomass values. In contrast, Rowbottom [START_REF] Rowbottom | Potential of Unmanned Aerial Vehicles (UAV) and Remote Sensing to Accurately Estimatepasture Biomass in Intensively Grazed Dairy Pastures; School of Earth and Environment[END_REF], in a comparison between the rising plate meter (RPM) height and UAS as monitoring tools for pasture biomass, reached a lower modeling quality with the UAS imagery than with the RPM height. It is worth noting that this last study only used spectral information from a consumer grade digital camera and that their results in terms of variance explanation were similar to the ones found for the spectral biomass model in this study (Figure 4, Model 3). Our approach, which used both sward height and spectral information, allowed us to take advantage of a top canopy parameter (sward height) and subcanopy parameters, reflecting the activity of the swards over most of their depths (spectral data).

Mapping Sward Height Differences and Biomass with UAS

Our approach highlights pasture biomass spatial heterogeneity, which cannot be easily captured using classical pasture field monitoring tools for research or grazing management purposes. Capturing this spatial heterogeneity is essential as it plays a crucial role in animal grazing behavior and the selection of feeding stations across a given area.

The mapping of sward height and sward height differences at such a fine spatial resolution scale are essential to accommodate the full potential of virtual fencing, and they could give new perspectives to more traditional grazing management, such as the quick response to the adaptation of grazing schedules, provided that adequate automated data acquisition and treatment applications and decision support tools are developed for farmers.

Discussion

Modeling Biomass of Pasture with UAS Imagery

The performance (adj. R 2 = 0.49) of the best performing UAS-based biomass model (Model 4) was similar to the state-of-art field pasture biomass non-destructive monitoring tools for swards-turf grasses such as the rising plate meter [START_REF] Harmoney | Determination of pasture biomass using four indirect methods[END_REF][START_REF] Laca | Comparison of weight estimate and rising-plate meter methods to measure herbage mass of a mountain meadow[END_REF]. There are still very few references for the application of UAS imagery products to pasture biomass modeling on such intensively managed pastures with rather low pasture heights. Lee et al. [START_REF] Lee | Mapping Herbage Biomass on a Hill Pasture using a Digital Camera with an Unmanned Aerial Vehicle System[END_REF] achieved a better model performance (R 2 = 0.77) than this study, but their work was performed on a multisite dataset in Korea, which had a higher diversity and more contrasting biomass values. In contrast, Rowbottom [START_REF] Rowbottom | Potential of Unmanned Aerial Vehicles (UAV) and Remote Sensing to Accurately Estimatepasture Biomass in Intensively Grazed Dairy Pastures; School of Earth and Environment[END_REF], in a comparison between the rising plate meter (RPM) height and UAS as monitoring tools for pasture biomass, reached a lower modeling quality with the UAS imagery than with the RPM height. It is worth noting that this last study only used spectral information from a consumer grade digital camera and that their results in terms of variance explanation were similar to the ones found for the spectral biomass model in this study (Figure 4, Model 3). Our approach, which used both sward height and spectral information, allowed us to take advantage of a top canopy parameter (sward height) and subcanopy parameters, reflecting the activity of the swards over most of their depths (spectral data).

Mapping Sward Height Differences and Biomass with UAS

Our approach highlights pasture biomass spatial heterogeneity, which cannot be easily captured using classical pasture field monitoring tools for research or grazing management purposes. Capturing this spatial heterogeneity is essential as it plays a crucial role in animal grazing behavior and the selection of feeding stations across a given area.

The mapping of sward height and sward height differences at such a fine spatial resolution scale are essential to accommodate the full potential of virtual fencing, and they could give new perspectives to more traditional grazing management, such as the quick response to the adaptation of grazing schedules, provided that adequate automated data acquisition and treatment applications and decision support tools are developed for farmers.

Operational Recommendation

Considering the targeted spatial resolution of the UAS products and the scale of the study site (2 × 2000 m 2 plots), the choice of an octocopter seemed to be appropriate and is recommended for potential users of the developed approach. The use of fixed-wing UAS would induce higher and faster flights, with a negative impact on the spatial resolution of the UAS imagery and the quality of the 3D reconstruction. Moreover, the relatively low shutter speed of multispectral sensors (1/500 for the multispectral sensors of the Parrot Sequoia) can induce blurry images if acquired at higher flight speeds (>12 m.s -1 at 50 m above ground level for the Parrot Sequoia).

Our study leveraged the height information provided by the combination of photogrammetric DSM and a LiDAR digital terrain model acquired on a regional scale by public administration. The availability of aerial LiDAR for an entire region/country is becoming quite common in Western Europe, but is still rare in other parts of the world. As an alternative, a UAS flight survey after an early mowing of studied pastures growing could produce a reference initial digital surface model, which could then be used to normalize the DSM acquired later. Moreover, the computation of a high resolution photogrammetric DSM requires a long computational time and more specific photogrammetric expertise. Although these elements could reduce the replication of our approach in diverse contexts, they can be circumvented by the use of UAS LiDAR. Bringing a LiDAR sensor onboard would allow the direct production of a height model while decreasing the computing time and error in terms of sward height estimation. The cost of such an UAS LiDAR approach is still prohibitive for many, but one can assume that prices will drop in the near future.

The regression models proposed in our study must be calibrated if used in different contexts (species, season, stocking rates, etc.). Nonetheless, they still represent a valuable improvement compared to the use of state-of-art field approaches regarding model performance and the production of high resolution maps.

Conclusions

Our results demonstrate the potential of UAS imagery to be used as a tool for precision grazing study applications. The best performing pasture biomass model combined UAS sward height information with UAS vegetation indices. The quality of the fitted models (sward height or biomass) using UAS imagery opens up a new area of research in precision grazing. The UAS approach of height and biomass monitoring was revealed to be a potential alternative to more conventional time-consuming field methods. In addition, the quality of the models based on UAS imagery, UAS spectral information, and 3D imagery can be used to derive sward height differences and biomass maps at a very fine spatial scale. Considering the versatility of the UAS, future research should integrate pasture monitoring at a very fine spatial scale using denser time series at a sub-daily (even hourly) scale. As the products of UAS imagery are georeferenced precisely, they can be compared with other types of spatial information to monitor the grazing of single individuals. This study also demonstrated that 3D information (sward height) derived from the UASs could be useful for pasture biomass modeling. Even if it requires additional photogrammetric skills, UAS 3D information is a direct product of the photogrammetric process used to obtain the spectral orthoimages. The use of UAS 3D data should be promoted in further studies, as this information has been used relatively less often in agricultural studies, while other research fields, such as forestry, have been evaluating its use quite extensively. The use of other modeling approaches, such as mechanist models, could further improve the reliability of predicted and mapped pasture biomass results.

In terms of grazing behavior research, our approach could be combined with information provided by motion capture collars, such as those described by Andriamandroso et al. [START_REF] Andriamandroso | Development of an open-source algorithm based on inertial measurement units (IMU) of a smartphone to detect cattle grass intake and ruminating behaviors[END_REF]. Regarding the very fine spatial scale of the UAS products, attention must be paid to the quality of the positioning associated with the information to which it is compared. The recent arrival on the market of low-cost centimetric precision GPS and inertial measurement units will soon allow for the development of new motion capture systems, which combined with UAS imagery, will provide new research opportunities in terms of precision grazing and the characterization of grazing animal behavior.

Figure 1 .

 1 Figure 1. Study site location in Belgium (Western Europe) and locations of the field measurements (biomass sampling and terrestrial LiDAR scans)). P1 and P2 were grazed for two consecutive days.

  Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 14 sward height models. A cross validation of Model 1 was performed using the repeated k-fold approach (k = 5; repetitions = 100).

Figure 2 .

 2 Figure 2. Sampling scheme for the UAS sward height model (SHM) validation with the LiDAR sward height model. The SHMs were compared within a circular area of interest (2 m < radius < 5 m), following a 0.5 m square grid.

Model 2 : 2 )Model 3 : 3 )Model 4 :

 22334 Pasture biomass = f (UAS height) (Pasture biomass = f (UAS Reflectance, UAS VI) (Pasture biomass = f (UAS height, UAS Reflectance, UAS VI).

Figure 2 .

 2 Figure 2. Sampling scheme for the UAS sward height model (SHM) validation with the LiDAR sward height model. The SHMs were compared within a circular area of interest (2 m < radius < 5 m), following a 0.5 m square grid.

2. 4 . 2 )Model 3 : 3 )Model 4 :

 42334 Modeling Biomass of Pasture with UAS Imagery The modeling of pasture biomass (before grazing) was investigated comparatively with different UAS products (3D and/or spectral): Model 2: Pasture biomass = f (UAS height) (Pasture biomass = f (UAS Reflectance, UAS VI) (Pasture biomass = f (UAS height, UAS Reflectance, UAS VI).

Figure 3 .

 3 Figure 3. The linear regression (Model 1) between the reference sward height data (LiDAR) and UAS sward height data had good agreement between the two methods. The dashed line represents the fitted model and the solid line represents the 1:1 line (i.e., x = y).

Figure 3 .

 3 Figure 3. The linear regression (Model 1) between the reference sward height data (LiDAR) and UAS sward height data had good agreement between the two methods. The dashed line represents the fitted model and the solid line represents the 1:1 line (i.e., x = y).

Figure 4 .

 4 Figure 4. Comparison of different approaches of pasture biomass modeling with UAS. The Model 2 (I) is a biomass model based on UAS height. Model 3 (II) and Model 4 (III) are multilinear models using UAS spectral data (Model 3) and a combination of UAS spectral and height data (Model 4). For Model 2, the dashed lines represent the fitted models. The solid lines represent the 1:1 lines (i.e., x = y).

Figure 4 .

 4 Figure 4. Comparison of different approaches of pasture biomass modeling with UAS. The Model 2 (I) is a biomass model based on UAS height. Model 3 (II) and Model 4 (III) are multilinear models using UAS spectral data (Model 3) and a combination of UAS spectral and height data (Model 4). For Model 2, the dashed lines represent the fitted models. The solid lines represent the 1:1 lines (i.e., x = y).

Figure 5 .

 5 Figure 5. Sward height before (I) and after (II) cattle grazing, sward height differences (III) (sward height 'before' minus sward height 'after') and associated Kernel density plot (IV). The majority of the pixel values were associated with positive values (intakes or trampling by cattle).

Figure 5 .

 5 Figure 5. Sward height before (I) and after (II) cattle grazing, sward height differences (III) (sward height 'before' minus sward height 'after') and associated Kernel density plot (IV). The majority of the pixel values were associated with positive values (intakes or trampling by cattle).

Figure 6

 6 Figure 6 allowed us to visualize the spatial pattern of the pasture biomass and its heterogeneity at a very fine scale. This map was produced through the use of Model 4 applied to the associated UAS variables (i.e., UAS sward height, GRVI, GNDVI, and NDRE). The northwestern part of the study site presented lower biomass values. Some abnormal values were shown on the empirical cumulative distribution function (ECDF) plot (Figure6II) that were well above the usual averages observed for such pastures or were even negative. Nevertheless, these values can be considered as noise when considering that the 5th and 95th percentiles of the predicted (and mapped) pasture biomass values were 0.04 kg/m 2 and 0.36 kg/m 2 , respectively.

Figure 6 .

 6 Figure 6. Mapping pasture biomass before grazing using predictors of the best model (Model 4) within pastures P1 and P2 (I) and the empirical cumulative distribution function (ECDF) plot of the UAS pasture biomass values (II).

Figure 6 .

 6 Figure 6. Mapping pasture biomass before grazing using predictors of the best model (Model 4) within pastures P1 and P2 (I) and the empirical cumulative distribution function (ECDF) plot of the UAS pasture biomass values (II).

Table 2 .

 2 Variables extracted from the UAS imagery.

	Name	Type	Ground Sampling Distance of the Layer (m)	Sensor
	NDVI	Vegetation index (VI)	0.05	Multispectral (Sequoia)
	NDRE	Vegetation index (VI)	0.05	Multispectral (Sequoia)
	GNDVI	Vegetation index (VI)	0.05	Multispectral (Sequoia)
	GRVI	Vegetation index (VI)	0.05	Multispectral (Sequoia)
	Red (R)	Reflectance	0.05	Multispectral (Sequoia)
	Green (G)	Reflectance	0.05	Multispectral (Sequoia)
	Near Infra-Red (NIR)	Reflectance	0.05	Multispectral (Sequoia)
	Red-Edge (RE)	Reflectance	0.05	Multispectral (Sequoia)
	Sward Height Model (SHM)	3D	0.025	RGB (Sony RX100)

Table 3 .

 3 Relative importance (% of R²), estimates of coefficients, and p-values of the fitted linear model with the selected variables computed for Model 4, which combines the UAS sward height and UAS spectral data.

		Rel. Importance (%)	Regression Coefficients	Pr (>|t|)
	UAS height (SHM)	45	1.1	0.00105	**
	GRVI	27	-0.1	0.00139	**
	GNDVI	17	6.5	0.01162	*
	NDRE	11	0.8	0.04675	*
	Intercept	/	-4.2	0.01725

*

3.3. Mapping Sward Height Differences and Biomass with UAS

Table 3 .

 3 Relative importance (% of R 2 ), estimates of coefficients, and p-values of the fitted linear model with the selected variables computed for Model 4, which combines the UAS sward height and UAS spectral data.

		Rel. Importance (%)	Regression Coefficients	Pr (>|t|)
	UAS height (SHM)	45	1.1	0.00105	**
	GRVI	27	-0.1	0.00139	**
	GNDVI	17	6.5	0.01162	*
	NDRE	11	0.8	0.04675	*
	Intercept	/	-4.2	0.01725	*

Funding: This research received no external funding.

Acknowledgments:

The authors would like to thank the 'drone' team who performed the UAS flight surveys (particularly Cédric Geerts and Samuel Quevauvillers).