Valentin Chardon
email: valentin.chardon@live-cnrs.unistra.fr

Laurent Schmitt

Hervé Piégay

Dimitri Lague

Use of terrestrial photosieving and airborne topographic LiDAR to assess bed grain size in large rivers: a study on the Rhine River

Keywords: Grain size, Large rivers, above-water and underwater conditions, Methodological study, Terrestrial photosieving, Airborne topographic LiDAR

Most grain size monitoring is still being conducted by manual sampling in the field, which is time consuming and has low spatial representation. Due to new remote sensing methods, some limitations have been partly overcome, but methodological progress is still needed for large rivers as well as in underwater conditions. In this paper, we tested the reliability of two methods along the Old Rhine River (France/Germany) to estimate the grain size distribution (GSD) in above-water conditions: (i) a low-cost terrestrial photosieving method based on an automatic procedure using Digital Grain Size (DGS) software and (ii) an airborne LiDAR topobathymetric survey. We also tested the ability of terrestrial photosieving to estimate the GSD in underwater conditions. Field pebble counts were performed to compare and calibrate both methods. The results showed that the automatic procedure of 1 2 terrestrial photosieving is a reliable method to estimate the GSD of sediment patches in both above-water and underwater conditions with clean substrates. Sensitivity analyses showed that environmental conditions, including solar lighting conditions and petrographic variability, significantly influence the GSD from the automatic procedure in above-water conditions. The presence of biofilm in underwater conditions significantly altered the GSD estimation using the automatic procedure, but the proposed manual procedure overcame this problem. The airborne LiDAR topographic survey is an accurate method to estimate the GSD of above-water bedforms and is able to generate grain size maps. The combination of terrestrial photosieving and airborne topographic LiDAR methods is adapted to assess the GSD along large rivers in entire sections that are several kilometers long.

Introduction

Quantification of the GSD of a riverbed surface is important for assessing sediment transport and for geomorphological and ecological studies, notably regarding river restoration (Graham, Rice, and Reid, 2005b;[START_REF] Mackenzie | Breaking from the Average: Why Large Grains Matter in Gravel Bed Streams[END_REF].

Many studies have been carried out during recent decades to characterize the GSD, but its assessment remains difficult because the grain size variability depends notably on the study scale (Church, Mclean, & Wolcott J.F, 1987;Graham et al., 2005b). GSD assessment of a sediment patch requires a large set of particles whose number depends on the shape of the distribution [START_REF] Rice | Sampling surficial fluvial gravels: the precision of size distribution percentile estimates[END_REF]. Most grain size monitoring was performed by using the classic "pebble count" manual procedure [START_REF] Wolman | From Manned to Unmanned Aircraft: Adapting Airborne Particle Size Mapping Methodologies to the Characteristics of SUAS and SfM[END_REF]. These different procedures are all time consuming and show low spatial representation, especially for large rivers [START_REF] Warrick | Cobble Cam: Grain-Size Measurements of Sand to Boulder from Digital Photographs and Autocorrelation Analyses[END_REF]).

2

From the beginning of the 1970s, terrestrial photosieving methods were carried out to estimate the GSD in above-water conditions to reduce field time, increase the sampling density and limit bed disturbance [START_REF] Butler | Automated Extraction of Grain-Size Data from Gravel Surfaces Using Digital Image Processing[END_REF]. Data processing requires manually digitizing the B-axis of each visible particle, increasing the time and lab effort [START_REF] Kellerhals | Sampling Procedures for Coarse Fluvial Sediments[END_REF][START_REF] Ibbeken | Photo-sieving: A Method for Grain-size Analysis of Coarse-grained, Unconsolidated Bedding Surfaces[END_REF].

Thereafter, semiautomatic or automatic procedures were developed to reduce the processing time. They were based either on oriented-object approaches [START_REF] Mcewan | Estimating the Size Composition of Sediment Surfaces through Image Analysis[END_REF][START_REF] Sime | Information on Grain Sizes in Gravel-Bed Rivers by Automated Image Analysis[END_REF]Graham, Reid, and Rice, 2005a;[START_REF] Lucas | Comparison of Automated Grain Sizing of Gravel Beds Using Digital Images to Standard Grid and Random-Walk Pebble Counts[END_REF][START_REF] Detert | Automatic Object Detection to Analyze the Geometry of Gravel Grains-a Free Stand-Alone Tool[END_REF][START_REF] Strom | Comparison of Automated Image-Based Grain Sizing to Standard Pebble-Count Methods[END_REF][START_REF] Sulaiman | Catena Application of Automated Grain Sizing Technique ( AGS ) for Bed Load Samples at Rasil River : A Case Study for Supply Limited Channel[END_REF] or on statistical approaches [START_REF] Buscombe | An Automated and Universal Method for Measuring Mean Grain Size[END_REF][START_REF] Buscombe | Transferable Wavelet Method for Grain-Size Distribution from Images of Sediment Surfaces and Thin Sections, and Other Natural Granular Patterns[END_REF][START_REF] Sime | Information on Grain Sizes in Gravel-Bed Rivers by Automated Image Analysis[END_REF][START_REF] Warrick | Cobble Cam: Grain-Size Measurements of Sand to Boulder from Digital Photographs and Autocorrelation Analyses[END_REF]. However, this kind of method is suitable only for studies where prompt GSD data are needed.

Over the last two decades, new methods of remote sensing have been developed to estimate the GSD at larger spatial scales in above-water conditions [START_REF] Dugdale | Aerial Photosieving of Exposed Gravel Bars for the Rapid Calibration of Airborne Grain Size Maps[END_REF]. Some studies tested the reliability of airborne digital imagery processing, such as airborne multispectral [START_REF] Rainey | Mapping Intertidal Estuarine Sediment Grain Size Distributions through Airborne Remote Sensing[END_REF][START_REF] Black | Mapping Sub-Pixel Fluvial Grain Sizes with Hyperspatial Imagery[END_REF] predict the GSD is dependent on the resolution of the images, which must be commensurate with the particle size [START_REF] Lejot | Utilisation de La Télédétection Pour La Caractérisation Des Corridors Fluviaux : Exemples d ' Applications et Enjeux Actuels Characterisation of Alluvial Plains by Remote Sensing : Cases Studies and Current Stakes[END_REF]. Three-dimensional data acquired by terrestrial LiDAR or UAV photogrammetry were also used to predict the GSD by using the bed surface roughness as a proxy [START_REF] Heritage | Terrestrial Laser Scanning of Grain Roughness in a Gravel-Bed River[END_REF][START_REF] Brasington | Modeling River Bed Morphology, Roughness, and Surface Sedimentology Using High Resolution Terrestrial Laser Scanning[END_REF][START_REF] Storz-Peretz | Morphotextural Characterization of Dryland Braided Channels[END_REF][START_REF] Woodget | Subaerial Gravel Size Measurement Using Topographic Data Derived from a UAV-SfM Approach[END_REF][START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF]. The first empirical relation between the roughness height metric (rh) and the median grain size was found by Gomez (1993) from the bed surface profile and vertical photographs. These methods have the advantage that they simultaneously provide various information, such as roughness maps, channel morphological dynamics, and grain size maps [START_REF] Marteau | Application of Structure-from-Motion Photogrammetry to River Restoration[END_REF][START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF]. However, the application of terrestrial LiDAR and UAV photogrammetry techniques is not applicable in large study reaches in underwater conditions.

During the last decade, airborne topographic LiDAR techniques have been increasingly employed, notably in fluvial geomorphological and hydrologic studies [START_REF] Charlton | Application of Airborne Lidar in River Environments: The River Coquet, Northumberland, UK[END_REF][START_REF] Demarchi | Hierarchical Object-Based Mapping of Riverscape Units and in-Stream Mesohabitats Using Lidar and VHR Imagery[END_REF][START_REF] Bizzi | LiDAR-Based Fluvial Remote Sensing to Assess 50-100-Year Human-Driven Channel Changes at a Regional Level: The Case of the Piedmont Region, Italy[END_REF][START_REF] Shaker | Automatic Land-Water Classification Using Multispectral Airborne LiDAR Data for near-Shore and 46 River Environments[END_REF]. Dataset collection is performed by using the nearinfrared wavelength [START_REF] Mandlburger | Topo-Bathymetric LiDAR for Monitoring River Morphodynamics and Instream Habitats-A Case Study at the Pielach River[END_REF]. The data collection is fast and may cover large areas. This method provides not only accurate point clouds in three dimensions but also intensity values that correspond to the amplitude of the return signal [START_REF] Kashani | A Review of LIDAR Radiometric Processing: From Ad Hoc Intensity Correction to Rigorous Radiometric Calibration[END_REF]. [START_REF] Kukko | Effect of Incidence Angle on Laser Scanner Intensity and Surface Data[END_REF] showed that this metric is impacted partly by surface roughness, suggesting a probable link between the intensity values and the bed grain size. Some studies have been carried out to qualitatively estimate the GSD in underwater conditions, but they are rare due to technical limitations in constraining field conditions.

Qualitative underwater substrate mapping was produced by using side-scan sonar imagery [START_REF] Buscombe | Compositional Signatures in Acoustic Backscatter Over Vegetated and Unvegetated Mixed Sand-Gravel Riverbeds[END_REF][START_REF] Buscombe | Automated Riverbed Sediment Classification Using Low-Cost Sidescan Sonar[END_REF][START_REF] Hamill | Alluvial Substrate Mapping by Automated Texture Segmentation of Recreational-Grade Side Scan Sonar Imagery[END_REF]. [START_REF] Rüther | Verifying a Photogrammetric Method to Quantify Grain Size Distribution of Developed Armor Layers[END_REF] used the terrestrial photosieving method in underwater conditions along the Surna River to estimate the bed GSD following an automatic procedure but without a control dataset and calibration with field GSD. This paper addresses the inference of the GSD in large rivers from remote sensing in both above-water and underwater conditions. Two methods are tested in abovewater conditions: (i) terrestrial photosieving, a low-cost technique to rapidly estimate the bed GSD, and (ii) airborne topographic LiDAR data analysis using a laser wavelength λ=1064 nm. For underwater conditions, we tested the terrestrial photosieving method. This paper follows three objectives: (i) compare the reliability between the terrestrial photosieving technique and the airborne LiDAR topometric survey in above-water conditions, (ii) test the reliability of the terrestrial photosieving technique to estimate the GSD in underwater conditions and (iii) determine the parameters that can influence the accuracy of both techniques in both environmental conditions.

Materials and methods

Study area

The Rhine is the third largest river in Europe with a drainage basin of approximately 185,000 km² and a length of 1250 km. The Upper Rhine runs from Basel to Bingenam-Rhein (300 km long) into the Upper Rhine graben, where the hydrological regime is nivo-glacial (Figure 1.a). The mean annual discharge is 1059 m 3 /s at Basel [START_REF] Uehlinger | The Rhine River Basin[END_REF]. After two engineering phases of correcting and regulating the flow (Arnaud, 2012;[START_REF] Schmitt | Managing Flood Risk: Innovative Approaches from Big Floodplain Rivers and Urban Streams[END_REF]

Data collection

Collection of digital images

Between winter and autumn 2017, terrestrial photosieving surveys were carried out on both above-water and underwater patches by using a 12-MP Olympus TG-4 waterproof camera. The camera was placed in a plan view by using a telescopic rod.

The recorded elevation ranged from 0.6 m to 1.0 m, and the pixel size ranged from 0.18 mm to 0.30 mm in above-water conditions and from 0.24 mm to 0.40 mm in underwater conditions. A square gabarit with an area of 1 m² placed on the bed surface and a rigid rule were used for scale in above-water and underwater conditions, respectively. All the digital images recorded in above-water conditions were characterized as clean substrate (n=10), while the digital images collected in underwater conditions were composed of both clean substrate (n=4) and substrate with biofilm (n=6).

Airborne topographic LiDAR survey

In March 2017, airborne LiDAR data were acquired by using a Teledyne Optech Titan DW sensor. The flight was carried out at an altitude of 400 m along a reach of 35 km between KP 174,00 and KP 209.00. The Titan DW sensor is a topobathymetric sensor operating a near-infrared (NIR) laser (λ=1064 nm) typical of airborne topographic LiDAR and a green laser (λ=532 nm) with shallow-water bathymetry mapping capability [START_REF] Fernandez-Diaz | Capability Assessment and Performance Metrics for the Titan Multispectral Mapping Lidar[END_REF] In this study, only the NIR laser data were used on above-water bars.

The divergence angle of the NIR laser is 0.3 mrad, resulting in a laser spot size of approximately 0.12 m at nadir. The range density points and range mean distance between points, 38 to 49 points/m 2 and 0.08 m to 0.09 m, respectively, were calculated on three-point clouds located on three above-water deposits. To estimate the range noise of the flight, we estimated the variability of the z-values and intensity values on three ground control areas (sections of roads) by assuming that the difference in the point values should be equal to 0 on planar control areas [START_REF] Lague | Accurate 3D Comparison of Complex Topography with Terrestrial Laser Scanner: Application to the Rangitikei Canyon (N-Z)[END_REF]. The uncertainties of the survey along the z-axis were less than 2 cm and equal to 2.10 cm for the intensity values (Table 1 Table 1: Precision of the LiDAR topographic survey along three ground control areas.

LiDAR metrics were computed for a kernel size equal to 0.4 m.

Collection of field pebble counts

The collection of field pebble counts was based on the random surface sampling of 100 particles for each sediment patch during February 2017 and October 2017 [START_REF] Wolman | From Manned to Unmanned Aircraft: Adapting Airborne Particle Size Mapping Methodologies to the Characteristics of SUAS and SfM[END_REF]. A sampling of 21 pebble counts was collected in both above-water (n=18) and underwater conditions (n=3) in a large range of sediment patches (Figure 2). The low number of field samplings in underwater conditions resulted from adverse survey conditions due to bed imbrication, high water depth and flow velocity and the locally important thickness of fine sediments (greater than 0.1 m). On above-water surfaces, each particle was collected at every 2*Dmax (Dmax=maximum B-Axis)

spacing along a decameter. In underwater conditions, the particles were sampled following the random step-toe procedure [START_REF] Bunte | Sampling Surface and Subsurface Particle-Size Distributions in Wadable Gravel-and Cobble-Bed Streams for Analyses in Sediment Transport , Hydraulics , and Streambed Monitoring[END_REF]. A caliper was used to measure the B-axis of each particle. The distribution was truncated at 3 psi according to the results of [START_REF] Rice | Sampling surficial fluvial gravels: the precision of size distribution percentile estimates[END_REF], who demonstrated an underestimation for these particles during pebble count sampling. Each pebble count area was georeferenced by using a global positioning system (centimetric planimetric precision). For calibration purposes, ten field pebble counts in above-water 10 conditions and three field pebble counts in underwater conditions were compared to the terrestrial photosieving results. On the other hand, 18 field pebble counts were used to calibrate the LiDAR metrics. The ranges of the 50 th and 84 th percentiles of the sampled sediment patches were between 4.75 psi and 7.02 psi and between 5.04 psi and 7.49 psi, respectively.

Data processing

Automated processing of digital images

Digital Grain Size (DGS) software developed by [START_REF] Buscombe | Transferable Wavelet Method for Grain-Size Distribution from Images of Sediment Surfaces and Thin Sections, and Other Natural Granular Patterns[END_REF] was used to automatically estimate the GSD of the apparent particle B-axis. This software, which is based on the wavelet statistical approach, provided the distribution in grid by number form. From the recommendations of [START_REF] Barnard | Field Test Comparison of an Autocorrelation Technique for Determining Grain Size Using a Digital 'beachball' Camera versus Traditional Methods[END_REF] and Buscombe (2010), a control dataset composed of twenty digital images was used to evaluate the accuracy of the software. The control dataset was performed by using the grid tool implemented in ImageJ software.

Manual processing of digital images

A manual method was also tested to process the digital images to compare with the automatic procedure. The B-axis of particles was manually digitalized following a grid in which nodes were spaced by one Dmax to avoid sampling the same particle twice.

This was made by the same operator to minimize the subjective bias during the Baxis choice. This procedure was performed by using the grid tool implemented on ImageJ software. Sampling areas were equal to 100 times the area of the largest particle following [START_REF] Petrie | Statistical Approach to Sediment Sampling Accuracy[END_REF] recommendation.

11 2.3.3 Processing the airborne topographic LiDAR survey

Computation of the Rh metric

First, all vegetation areas were excluded from the point cloud to avoid computation errors by empirically filtering the intensity values because it can be a more accurate method than geometric classification approaches [START_REF] Reymann | Improving LiDAR Point Cloud Classification Using Intensities and Multiple Echoes To Cite This Version : HAL Id : Hal-01182604 Improving LiDAR Point Cloud Classification Using Intensities and Multiple Echoes[END_REF]. The Rh metric was computed by using the roughness tool implemented in CloudCompare v.2.9.1 software (EDF R&D). For each point, the computation of the Rh metric is

given by the distance between the point elevation and the ordinary least-squares best fitting plane computed on its nearest neighbors (Woodget, 2015;[START_REF] Woodget | Subaerial Gravel Size Measurement Using Topographic Data Derived from a UAV-SfM Approach[END_REF][START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF]. The nearest neighbors are selected within a specified kernel size from the selected point. Third, according to the methodology of [START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF], the point cloud was clipped by a buffer with a radius of five meters, and the center corresponded to the location of the field pebble counts to compare the Rh metric distribution with the field GSD. This step was performed by using ArcGIS software.

Computation of the standard deviation of the intensity values

The intensity values of the return pulses of the NIR wavelength are influenced by several parameters, including internal instrument amplification, environmental effects (wetness and atmospheric transmittance), data acquisition geometry (range and incidence angle) and surface characteristics (reflectance and roughness) (Li et al., 2014;[START_REF] Kashani | A Review of LIDAR Radiometric Processing: From Ad Hoc Intensity Correction to Rigorous Radiometric Calibration[END_REF]. In our case study, instrumental and environmental parameters had a limited impact on intensity values because the same instrument was used, and environmental conditions, including the flight elevation, were relatively constant during flight. Considering the data acquisition geometry, [START_REF] Kukko | Effect of Incidence Angle on Laser Scanner Intensity and Surface Data[END_REF] and [START_REF] Kaasalainen | Analysis of Incidence Angle and Distance Effects on Terrestrial Laser Scanner Intensity: Search for Correction Methods[END_REF] showed that intensity 12 values are significantly sensitive to the incidence angle when the latter is greater than 20°. In our case, the incidence angle did not exceed 16°, and we assumed that the data acquisition geometry had no influence. Hence, we tested the standard deviation of the raw values of intensity as a proxy for the surface GSD. The first step consisted of calculating the standard deviation of the raw intensity values for each point of the point cloud for a given area by computing a buffer with a specific radius. We performed this step by using ArcGIS software. The second step was identical to the third step of the computation of the Rh metric (see section 2.3.4).

Sensitivity analyses

Effect of environmental conditions on GSD estimation by the DGS software

Sensitivity analyses were performed on three parameters:

1. The influence of the textural variation of the digital image induced by petrographic variability. This study was performed by using a median filter that homogenized the textural differences within individual particles and reduced the oversegmentation [START_REF] Strom | Comparison of Automated Image-Based Grain Sizing to Standard Pebble-Count Methods[END_REF]. The kernel size of the median filter was increased progressively to 1%, 2.5%, and 5% for underwater conditions (n=6) and increased progressively to 1%, 2.5%, 5%, 10% and 15% of the maximal B-axis of the largest particle observed on each digital image for above-water conditions with controlled solar lighting conditions (n=10). The GSDs estimated by the software were compared to the control dataset.

2. The light conditions during digital image acquisition. We recorded two photographs for each sampling patch (n=10) with and without sun exposure by using an umbrella on the above-water clean substrate. The time recording was extracted from the raw data digital image to calculate the influence of the solar angle on the GSD estimation.

13

3. The quality of the substrate on the DGS error estimations as clean substrate

or substrate with the presence of biofilm.

Selection of minimal sampling area in the DGS software

To evaluate the effect of the minimal sampling area on the estimation fluctuations of the DGS software, for each sediment patch in above-water conditions (n=10), we increased the minimal sampling area progressively from 10 to 100 times the area of the largest particle by steps of 10 by using a montage of digital images performed with GIMP software. We calculated a mobile slope coefficient to determine from which sampling area the fluctuation estimations were low as follows:

Slope coefficient =¿ ¿ 2.4.3 Airborne topographic LiDAR
The accuracy of the GSD predicted by the Rh metric can be influenced by the kernel size computation of the metric [START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF][START_REF] Woodget | Subaerial Gravel Size Measurement Using Topographic Data Derived from a UAV-SfM Approach[END_REF]. A sensitivity analysis was performed by computing the Rh metric for kernel sizes ranging between 0.4 m and 2.0 m, by steps of 0.1 m from 0.4 m to 1.0 m, and by steps of 0.2 m from 1.0 m to 2.0 m. This procedure was performed using CloudCompare v.2.9.1 software. The same procedure was performed for the standard deviation of the intensity values using ArcMap v10.3 software by progressively increasing the size of the radius of the computed buffer.

We computed linear regressions between the two proxies with the field GSD for each kernel size to define the best relationship between the two variables and consequently to identify the best kernel size. 2). Error estimations decrease as the percentiles increase for these two types of ground digital images. The NRMSE values were 88% and 32% for the 5 th percentile and the 95 th percentile, respectively.

15 the median filter reduced the NRMSE value of the estimated percentiles by the DGS software up to a certain threshold. Figure 6.a shows that 5% of the B-axis of the largest particle was the best compromise to minimize the estimated errors for all percentiles. The application of this filter on the raw-clean substrate reduced the estimated error percentiles by twice to four times (Table 2) and significantly reduced the over-segmentation by the software (Figure 3). The NRMSE values decreased from 83% to 25%, from 72% to 29% and from 51% to 18% for the 16 th , 50 th and 84 th percentiles, respectively.

The use of a median filter increased the data processing time. As shown in Figure 6.b, the computation of the median filter increased with increasing kernel size following a power function. For example, the computation time of a digital image with a size equal to 2748*2808 pixels was processed from 7 to 95 seconds for a kernel size equal to 1% and 5% of the B-axis of the largest particle.

Effect of sampling area

For all percentiles, a significant negative exponential relationship (95% confidence level) was found between the slope coefficient and the dimensionless sampling area of Graham et al. (2005b) was transformed to grid by number distribution. (a) Mean irreducible error calculated on three rivers.

Accuracy of the GSD estimation in above-water conditions by using airborne topographic LiDAR data

Rh metric and standard deviation of intensity values as proxies of the GSD

The percentiles calculated using the distribution of the standard deviation of the intensity values were higher than the threshold uncertainties for all percentiles (Figure 9.a). For the Rh metric, only the percentiles equal to or greater than the 25 th were greater than the threshold uncertainties.

23 Significant power law correlations at the 95% confidence level were found between the percentiles of the Rh metric distribution and the GSD of the field pebble counts in a range of grain sizes equal to or higher than the 25 th percentile (Table 4). The strongest correlation was found for the 84 th percentile with an R 2 coefficient equal to 0.85 (p-value4.9e-08). Significant power law correlations at the 95% confidence level were also found between the percentiles of the distribution of the standard deviation of the intensity values and the GSD of the field pebble counts (Table 4).

The strongest correlation is also found for the 84 th percentile with an R 2 coefficient equal to 0.85 (p-value9.6e-07) (Table 4). Figure 10.a-b shows that the relationship between the predicted percentiles by both metrics and the estimated percentiles by the field pebble counts are close to the 1:1 line. The mean values of the RMSE and MAE of the predicted GSD were equal to 0.34 psi and 0.29 psi for the Rh metric for the 16 th and 95 th percentile, respectively, while for the standard deviation of the intensity values, the mean values of the RMSE and MAE were equal to 0.38 psi and 0.30 psi for the 5 th and 95 th percentile, respectively (Figure 10.a-b). coefficient from 0.72 to 0.62 was observed for a kernel size increasing from 0.1 to 2.0 m. The highest R 2 was found for a kernel size equal to 0.4 m (grain size ranging between the 25th and 95 th percentiles). For the percentiles ranging between the 5 th and 95 th percentiles, R 2 decreases from 0.62 to 0.58 when the kernel size increases The differences in the predictions between the two proxies are given in (c). The blue and red dashed lines correspond to the mean and +/-1 standard deviation of the differences between the two metrics, respectively. The results showed that the accuracy of the GSD estimated by DGS software was dependent on the type of substrate. Low values of the NRMSE were obtained for clean substrate, while the NRMSE values exceeded 40% for all percentiles for substrate with the presence of biofilm (Table 2). For both types of substrate, the NRMSE decreased when the percentiles increased. Low values of the NRMSE for clean substrate are explained by low particle intravariability and high particle intervariability, whereas high NRMSE values for substrate with the presence of biofilm and/or presence of deposition are due to both low particle intravariability and intervariability of gray values (Figure 5). Considering sensitivity analyses, the best kernel size for digital image computation with the presence of biofilm to significantly reduce error estimations by the DGS software is equal to 2.5% of Dmax (Figure 14, Table 2). In this case, the median filter increases the intervariability between the particles (Figure 5). However, the NRMSE values are still too high to use DGS software to estimate the GSD for substrate with biofilms (Table 2). 3, Table 4). This can be explained by the fact that photosieving estimates grain size in the plan view, which is less affected by sediment patch imbrication than three-dimensional views [START_REF] Carbonneau | Automated Grain Size Measurements from Airborne Remote Sensing for Long Profile Measurements of Fluvial Grain Sizes[END_REF][START_REF] Wolman | From Manned to Unmanned Aircraft: Adapting Airborne Particle Size Mapping Methodologies to the Characteristics of SUAS and SfM[END_REF].

The automatic procedure of the digital images systematically underestimated the GSD of the field pebble counts, as also reported by [START_REF] Sime | Information on Grain Sizes in Gravel-Bed Rivers by Automated Image Analysis[END_REF][START_REF] Strom | Comparison of Automated Image-Based Grain Sizing to Standard Pebble-Count Methods[END_REF] and [START_REF] Warrick | Cobble Cam: Grain-Size Measurements of Sand to Boulder from Digital Photographs and Autocorrelation Analyses[END_REF], but it can be corrected by simple linear regressions (Figure 8). This underestimation by photosieving is explained by burial, overlapping, foreshortening and the C-axis not being oriented perpendicularly to the plane of some particles [START_REF] Graham | Maximizing the Accuracy of Image-Based Surface Sediment Sampling Techniques[END_REF]. The irreducible errors found in this study are in the same range as those in other studies (Table 3).

Comparison between the two proxies extracted from the airborne topographic LiDAR survey shows that the standard deviation of intensity seems to predict the GSD with higher accuracy than does the Rh metric and allows an estimation of the low percentiles of the distribution, contrary to the Rh metric. We explain this by the fact that the intensity values are less affected by recording uncertainties (Table 1). The comparison of our results with previous studies shows a large variability in the empirical relationships between the Rh metric (D84) and GSD (D84) for the same range of grain sizes (Figure 12). Our empirical relation is located between the relations of [START_REF] Woodget | Subaerial Gravel Size Measurement Using Topographic Data Derived from a UAV-SfM Approach[END_REF] and [START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF]. [START_REF] Pearson | Can High Resolution 3D Topographic Surveys Provide Reliable Grain Size Estimates in Gravel Bed Rivers?[END_REF] showed that bed roughness is influenced by several parameters, including the shape of particles, packing and bed imbrication of particles, and distribution shape and survey errors, and can explain differences in roughness values for the same range of grain sizes. [START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF] reported that their sediment patches were composed of well-rounded and subspherical granitic and metamorphic gravels and cobbles. In our case, particles were mostly flat and imbricated, generating lower roughness values for a given GSD. [START_REF] Woodget | Subaerial Gravel Size Measurement Using Topographic Data Derived from a UAV-SfM Approach[END_REF] also reported that their sediment patches were mostly composed of packed and imbricated particles. As intensity values are dependent on the bed surface roughness, it is very likely that various relationships could also be found between the standard deviation of the intensity values and the bed grain size for the same range of GSD. The approach of [START_REF] Pearson | Can High Resolution 3D Topographic Surveys Provide Reliable Grain Size Estimates in Gravel Bed Rivers?[END_REF] in laboratory conditions with the objective of assessing the influence of the sedimentary characteristics on the roughness -grain size relationship could also be performed in the field to try to build an abacus and derive from it several universal calibration laws usable for the two methods based on LiDAR data. It should be noted that this technique is limited only to sections of rivers where the roughness is relatively high and higher than the georeferenced error of the survey along the z-dimension.

Both methods have their own advantages: terrestrial photosieving allows the GSD to be followed at a high temporal frequency with rapid grain size information and is cheaper, whereas airborne topographic LiDAR can predict the GSD in large above-34 water areas and produce grain size maps, which is crucial to surveying the grain size evolution, notably for restorations along large rivers (Figure 13).

Toward an optimization of the two methods

Our study reveals that the accuracy of the automatic procedure of digital images to estimate the GSD is dependent on both the environmental conditions and the sampling area. The higher NRMSE values in our study compared to those published by [START_REF] Buscombe | Transferable Wavelet Method for Grain-Size Distribution from Images of Sediment Surfaces and Thin Sections, and Other Natural Granular Patterns[END_REF] may be explained by the petrographic variability of the particles, which induces oversegmentation due to high intravariability of the gray values of pixels (Figure 5). However, the use of a median filter can significantly reduce the oversegmentation with a decrease in the NRMSE value, close to the results of [START_REF] Buscombe | Transferable Wavelet Method for Grain-Size Distribution from Images of Sediment Surfaces and Thin Sections, and Other Natural Granular Patterns[END_REF]. Moreover, our study shows that DGS estimations were influenced by solar lighting in the majority of the cases with an underestimation using sunny photographs by comparison with photographs with shadow conditions. Changes in underestimation values are explained by the changes in the solar angle during a day.

These observations are in accordance with [START_REF] Warrick | Cobble Cam: Grain-Size Measurements of Sand to Boulder from Digital Photographs and Autocorrelation Analyses[END_REF] and may be explained by the fact that bed roughness generated by the coarser particles induced large shadow areas in low sun angles [START_REF] Warrick | Cobble Cam: Grain-Size Measurements of Sand to Boulder from Digital Photographs and Autocorrelation Analyses[END_REF] and explain why low percentiles were more impacted by solar conditions than were coarse percentiles.

Thus, the variation in solar lighting conditions can alter the calibration step between the GSD estimated by digital images and the field pebble counts. In addition, our study reveals that sensitivity analyses must be conducted to define the minimal sampling area from which the fluctuation of the DGS software becomes low to minimize variability during the calibration step with the field pebble counts. For our sediment patches, our study reveals that fluctuations were significantly reduced for a sampling area 80 times larger than the largest particle area. This sampling area is 35 dependent on the sorting index of the sampling sediment patches (Figure 7), but other parameters of the distribution shape (e.g., skewness, bimodality) can also influence the accuracy for a given surface sampling [START_REF] Rice | Sampling surficial fluvial gravels: the precision of size distribution percentile estimates[END_REF].

For the airborne topographic LiDAR, our study highlights that the accuracy of the Rh metric to predict the GSD can be dependent on the kernel size computation, contrary to the standard deviation of the intensity values, which requires no presensitivity analysis and significantly reduces the processing time (Figure 11).

The best kernel size to compute the Rh metric is four times larger than Dmax, which is twice the values found by [START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF]. We hypothesize that this difference is the consequence of a lower point density in our study (38 to 49/m 2 ) compared to that of [START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF]36 to 1052/m 2 ).

Supplementary research is needed to evaluate the impact of the density of the point cloud on the reliability of the bed GSD estimation by the intensity metric. [START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF] showed that 10 pts/m 2 is the minimal threshold to estimate the bed GSD by the Rh metric with accuracy for their range of sediment sizes. Research must also be conducted to evaluate the influence of the spot size of the beam, which depends both on the flight elevation and on the open angle of the captor.

Underwater conditions: accuracy of digital images

Our results show that the automatic procedure is a reliable method to estimate the GSD for clean substrate. Contrary to above-water conditions, there is no need to apply a median filter to reduce estimation errors by the DGS software. This is explained by the absorption of solar light by the water column, which significantly reduces the variability of the gray values on the digital images, thus minimizing the oversegmentation by the DGS software. However, this method does not allow substrate with biofilms to be processed even with the application of a median filter for percentiles equal to and greater than the 50 th percentile. An alternative approach must be conducted to process this kind of data, as this situation is very frequent for rivers with weak bedload dynamics and/or with high wash loads. The manual approach proposed in this paper is a reliable method to estimate the GSD of substrate with biofilms and has the advantages of requiring no preprocessing step and is independent of environmental conditions (Figure 15). This approach can also be applied in above-water conditions if photosieving was conducted in limiting environmental conditions (e.g., solar light for calibration purposes with field pebble counts, local growth of small pioneer vegetation). Unfortunately, it presents a nonnegligible operator-introduced bias during the manual B-axis determination.

Methodological recommendations

To perform data calibration from terrestrial photosieving with field pebble counts in above-water conditions, we recommend the following:

(i) Collect the dataset during cloudy days as proposed by [START_REF] Warrick | Cobble Cam: Grain-Size Measurements of Sand to Boulder from Digital Photographs and Autocorrelation Analyses[END_REF][START_REF] Buscombe | An Automated and Universal Method for Measuring Mean Grain Size[END_REF] and when the vegetation is not yet developed (Graham et al., 2005b).

(ii) Use a median filter in case of large variability in particle petrography to reduce the oversegmentation by the DGS software. The best kernel size for the median filter computation found in this study can vary as a function of the petrographic characteristics of the sediment patches studied; thus, new sensitivity analyses must be conducted for studies performed on other rivers.

(iii) Perform sensitivity analyses to evaluate the minimal sampling area from which the estimation fluctuations of the DGS software are judged to be low in relation to the characteristics of sediment patches studied.

37

We used between ten and twenty field data calibrations for the two methods used in this study deployed along several kilometers of river reach. However, the number of these survey calibrations may be higher for studies covering several tens of kilometers because the heterogeneity of the GSD may be relatively important due to various sedimentary contexts (e.g., influence of tributaries, bedload inputs, contacts with lateral hillslope), which can generate various sediment patch characteristics [START_REF] Vázquez-Tarrío | Using UAS Optical Imagery and SfM Photogrammetry to Characterize the Surface Grain Size of Gravel Bars in a Braided River (Vénéon River, French Alps)[END_REF].

Conclusion

This study highlights the ability of terrestrial photosieving and, for the first time, the capacity of airborne topographic LiDAR surveys based on the computation of two metrics, the Rh metric and the standard deviation of the intensity values, to estimate the GSD of above-water deposits. The errors of both methods are similar with a range of MAE values from 0.23 to 0.31 psi. Airborne topographic LiDAR allows the production of maps of the GSD, while rapid GSD estimations can be produced only by terrestrial photosieving. The results also show that the accuracy of automatic processing of terrestrial photosieving depends on the environmental conditions and areal sampling. For both proxies derived from airborne topographic LiDAR, the accuracy of the predicted GSD for the Rh metric is dependent on the size of the kernel computation, contrary to the standard deviation of the intensity values. The latter allows the estimation of all grain size percentiles, whereas the Rh metric allows the estimation of only the grain size percentiles coarser than the 25th due to the larger uncertainties of recording. The automatic procedure of terrestrial photosieving is also a reliable method to estimate the GSD in underwater conditions for clean substrates, but it cannot be adapted in the presence of biofilm. The manual procedure can be a relevant alternative to overcome these constraints. Thus, combining terrestrial photosieving with airborne topographic LiDAR appears a relevant compromise for surveying the GSD in large rivers, along reaches of several kilometers and over the whole channel cross section with rapid GSD data in underwater conditions and continuous GSD data in above-water conditions. In the future, one important challenge is to develop universal laws according to sediment patch characteristics for both terrestrial photosieving and airborne topographic LiDAR surveys.

  or color photography, by performing textural analyses to estimate the GSD over scales from several hundred meters to several dozen kilometers (The accuracy of these methods to

Figure 1 :

 1 Figure 1: (a) The Rhine Basin and the location of the study reach, (b) the Old Rhine River and the GCA, (c) the range of digital images collected in this study, and (d) and (e) maps of the intensity values and the Rh roughness metric of an above-water deposit along the study reach. In (a), the black rectangle indicates the location of the Old Rhine River. In (b), the black arrows, the black cross and the gray circles

Figure 2 .

 2 Figure 2. (a) Framework of the study and maximum daily flow at the Rheinweiler gauging station during 2017 and (b) the GSD of the field pebble counts (truncated at 3 psi). * indicates the collection period of the field pebble counts used for the

1

 1 Figure3shows that the percentiles estimated by the DGS software in clean abovewater conditions with clean substrate and controlled solar lighting conditions were systematically underestimated compared with the percentile values of the control dataset. For the 50 th percentile, the values of the NRMSE were equal to 72% for the digital images with clean substrate in above-water conditions (Table2). Error

Figure 3 :

 3 Figure 3: Estimated percentiles by the DGS software as a function of the percentile values of the control dataset. The dashed lines correspond to the 1:1 line. The digital images were taken with controlled solar lighting conditions in above-water conditions.

Figure 4

 4 Figure4compares the percentiles estimated by the DGS software for the same digital image in sunny and shadow conditions. An overprediction below 10% was observed before 10:00 am for sunny conditions by comparison with shadow conditions. In contrast, underprediction was observed in sunny conditions with a peak of 20% for 2:00 pm. In any light condition, coarse percentiles are less impacted than small percentiles.

Figure 4 :

 4 Figure 4: Impact of the solar lighting conditions on the GSD estimated by DGS software.

Figure 5 :

 5 Figure 5: Effect of applying a median filter on gray values of digital images in different conditions: (a) and (d) correspond to a digital image of clean substrate in controlled solar lighting conditions in above-water conditions with and without the application of a median filter, respectively; (b) and (e) correspond to a digital image of clean substrate in underwater conditions with and without the application of a median filter, respectively; (c) and (f) correspond to a digital image of substrate with the presence of biofilm in underwater conditions with and without the application of a median filter, respectively. The gray line in (g), (h) and (i) corresponds to the gray values of pixels

Figure 6 :

 6 Figure 6: (a) Effect of applying the median filter to the model deviation for each percentile as a function of the kernel size. Percentage values correspond to the kernel size of the median filter divided by the B-axis of the largest particle, multiplied by 100.

(Figure 7 :

 7 Figure 7: (a) Evolution of the slope coefficient as a function of the dimensionless sampling area and the sorting index, where black lines correspond to the relationship, and (b) boxplots of the slope values according to the dimensionless sampling area. Intervals correspond to the 25 th and 75 th percentiles, and the lower and upper whiskers correspond to 1.5 times the interval quartile range. Black lines correspond to the median values, and the red points correspond to the mean values.

Figure 8 :

 8 Figure 8: Percentiles estimated by field pebble counts according to percentiles estimated by the automatic procedure of digital images (a) before the bias correction and (b) after the bias correction. RMSE and MAE values are in psi. The distributions were truncated at 3 psi.

Figure 9 :

 9 Figure 9: Values of percentiles estimated by (a) the standard deviation of the intensity values and (b) the Rh metric extracted from the airborne topographic LiDAR survey. The black line corresponds to the threshold uncertainties of the two metrics estimated on the three control areas (see section 2.2.2).

Figure 10 :

 10 Figure 10: Predicted percentiles by (a) the Rh metric and (b) the standard deviation of the intensity values for different percentiles.

  from 0.1 to 2.0 m. No clear threshold of kernel size appears for the standard deviation of the intensity, contrary to the Rh metric.

Figure 11 :

 11 Figure 11: Impact of the kernel size used for the computation of the two proxies for (a) the Rh metric and (b) the standard deviation of the intensity values on the statistical relationships between the distribution of the metrics and the GSD of the field pebble count. Black points correspond to the mean values of R 2 from the 25 th to the 95 th percentiles; the gray points correspond to the mean values of R 2 for the 5 th to the 95 th percentiles. Error bars correspond to +/-one standard deviation.

Figure 12 :

 12 Figure 12: Relation between the 84 th percentile of the field pebble-count distribution according to the 84 th percentile of the distribution of the Rh metric compared with the literature results. The calibration of Woodget et al., (2017) was based on a grid count on ground digital images (0.4 m x 0.4 m), whereas Vázquez-Tarrío et al., (2017) used field pebble counts. Values of correspond to the slope of the linear relations found between the two variables.

Figure 13 :

 13 Figure 13: Predicted maps of the median grain size on an above-water bar by using (a) the standard deviation of the intensity values and (b) the Rh metric as a proxy.

  photosieving

Figure 14 : 1

 141 Figure 14: Effect of the application of the median filter on the NRMSE values for each percentile according to the kernel size computation. Percentages correspond to the kernel size of the median filter divided by the B-axis of the largest particle multiplied by 100.

  ).

	Standard deviation of z	Standard deviation of Rh		Standard deviation of
			Standard deviation of
	values deviation	metric deviation (cm)		intensity metric deviation
			intensity values deviation
	(cm)
	1.00	0.18	2.10	6.51
	2.00	0.17	1.90	3.65
	0.90	0.15	2.10	1.50
		9

Table 2 :

 2

		Above-water condition		Underwater condition
	Decile	Clean substrate with	Clean substrate with	Clean substrate	Biofilm	Biofilm with median
			16

Values of the normalized root mean square error (NRMSE) in percentage for each type of data. The RMSE values were normalized by the mean particle size of the corresponding percentile of the control dataset and multiplied by 100. 3.1.2 Effects of environmental conditions 3.1.2.1 Effects of solar lighting conditions

Table 3 :

 3 Table 3). Comparison of the performance of the procedures used in this study with those used in other published studies. All values are in psi. The irreducible error (e)

				Form	Lighting	Range of grain		Irreducible
			Approach				Truncations
	Study			distribution	conditions	size (psi)		error (e) (psi)
	This study		Statistical	Pebble	Controlled	3.00 -8.13	3 psi	0.29
			approach	counts
	This study		Manual	Pebble	Not controlled	3.00 -8.13	3 psi	0.29
			approach	counts
	Sime	&
			Object-oriented	Pebble
	Ferguson,				Not specified	Not specified	3 psi	0.27
			approach	counts
	(2003)
	Graham et al., (2005.b)	Object-oriented approach	Pebble counts	Controlled	0.10 -7.64	3 psi	0.18	(a)
	Strom et al.,	Object-oriented	Pebble	Not specified	2.32 -7.60	2.32 psi	0.26
	(2010)		approach	counts

Table 4 :

 4 Statistical relationships between the percentiles of the distributions of the Rh metric and the standard deviation of the intensity values and the GSD of the field pebble counts. The RMSE and MAE metrics are in psi. NE corresponds to no estimation.

			Rh metric				Standard deviation of intensity values
			(n=18)						(n=18)
	Percentile	Equation	R 2	p-value	RMSE	MAE	Equation	R 2	p-value	RMSE	MAE
	D5	NE	NE	NE	NE	NE	y2.528x 0.154	0.32	0.013	0.40	0.31
	D10	NE	NE	NE	NE	NE	y2.41x 0.180	0.41	0.0048	0.39	0.31
	D16	NE	NE	NE	NE	NE	y12.138x 0.219	0.61	0.00035	0.32	0.28
	D25	y30.280x 0.345	0.62	9.4e-05	0.36	0.30	y1.800x 0.264	0.67	0.00015	0.34	0.30
	D50	y23.994x 0.319	0.69	2.2e-05	0.38	0.31	y1.450x 0.318	0.78	3.2e-05	0.32	0.26
	D75	y23.654x 0.341	0.79	8.8e-07	0.34	0.28	y1.110x 0.387	0.82	2.2e-06	0.33	0.25
	D84	y29.015x 0.410	0.85	4.9e-08	0.29	0.26	y1.103x 0.386	0.85	9.6e-07	0.31	0.24
	D90	y28.590x 0.422	0.81	3.5e-07	0.33	0.32	y1.096x 0.386	0.83	1.9e-06	0.33	0.23
	D95	y24.781x 0.398	0.84	8.0e-08	0.29	0.26	y1.35x 0.338	0.76	8.9e-06	0.36	0.28

Acknowledgements

This research was supported by the company Electricité de France (EDF) within the research collaboration 'Redynamization of the Old Rhine -years 2014-2018' (EDF 5910132058). The Teledyne Optech Titan DW lidar sensor used in this study is operated by the Nantes-Rennes Topo-Bathymetric lidar platform (University of Rennes 1/University of Nantes). The authors thank Patrice Carbonneau for helpful recommendations and Pierre-Alexis Herrault for the computation of metrics characterizing the point clouds.

To perform data calibration from the two metrics extracted from the airborne topographic LiDAR survey with the field pebble counts in above-water areas, we recommend the following:

(i) Prioritize the dataset collection in winter due to the absence of vegetation to limit estimation errors during the computation of the two metrics and evaluate the GSD on the largest possible area. Preferentially use the standard deviation of the intensity values rather than the Rh metric because the accuracy prediction of the metric is weakly dependent on the kernel size computation of the metric and reduces processing time.

(ii) In the case when intensity values are not available, the Rh metric can be a useful proxy, but sensitivity analyses must be performed to find the best size of the kernel computation of the metric.

The relationships found in this study have been performed for a maximum beam size equal to 12 cm and a minimal density of points greater than 30 pts/m 2 . Supplementary research should be conducted to evaluate the influence of these parameters on the relation between the two proxies tested in this study and the bed GSD.

For digital image collection in underwater conditions, we recommend the following:

(i) Perform a survey once the water turbidity is low and biofilm is not yet developed as during low water temperature or just after critical flood events.

(ii) Collect the data by using the video option of the camera rather than the photographs option to limit blurred images due to the movement of the boat, which does not allow the captor to make the adjustment and requires less time in the field.


David Dellong

Frauke Klingelhoefer

Anke Dannowski

Heidrun Kopp

Shane Murphy

David Graindorge

Lucia Margheriti

Milena Moretti

Giovanni Barreca

Luciano Scarfì

Alina Polonia

Marc-Andre Gutscher

Reply to Comment by A. Argnani on "Geometry of the Deep Calabrian Subduction from Wide-Angle Seismic Data and 3-D Gravity Modeling"

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

Introduction

First of all, we would like to thank Andrea Argnani for the interest he has shown in our work. Andrea Argnani in his comment on Dellong et al., 2020 (Geometry of the deep Calabrian subduction (Central Mediterranean Sea) from wide-angle seismic data and 3-D gravity modeling), proposes an alternate interpretation of the final velocity models presented by [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF]Dellong et al., 2020. Additionally, he claims that a part of the literature was not reported properly. While we basically agree with some of the interpretations presented in his comment (activity of the northern Malta Escarpment, presence of subducting crust at the base of the Calabrian crustal block), we fundamentally disagree with him on the presence and the activity of the lithospheric tear fault in the Ionian Sea south of Mt. Etna. We will respond to his points in the order that he raised them.

Active tectonics and Tear faults at the Malta Escarpment.

Tectonic activity of the Malta Escarpment (ME)

The activity of the ME has been discussed previously [START_REF] Gallais | Propagation of a lithospheric tear fault (STEP) through the western boundary of the Calabrian accretionary wedge offshore eastern Sicily (Southern Italy)[END_REF][START_REF] Gutscher | Tectonic expression of an active slab tear from high-resolution seismic and bathymetric data offshore Sicily (Ionian Sea)[END_REF][START_REF] Gutscher | Active tectonics of the Calabrian subduction revealed by new multi-beam bathymetric data and high resolution seismic profiles in the Ionian Sea (Central Mediterranean)[END_REF][START_REF] Polonia | The Ionian and 205 Alfeo-Etna fault zones: New segments of an evolving plate boundary in the central Mediterranean Sea?[END_REF][START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF]. We agree with Argnani that the ME is a currently active structure in its northern portion (North of Siracusa) showing primarily normal faulting (extension to transtension). However, the central and southern part of the ME (south of Siracusa) does not exhibit significant crustal earthquakes nor seismic images of active faulting. This aspect of Argnani's comment seems much more directed towards earlier work [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF], which clearly explains why the Alfeo Fault System is preferred as the current-day location of the STEP (Subduction Transform Edge Propagator). [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF] specify that the Malta escarpment affects the crust but represents probably an inherited structure from the Early Mesozoic and previous studies show that activity is mostly normal [START_REF] Torelli | Plio-Quaternary tectonic evolution and structure of the Catania foredeep, the northern Hyblean Plateau and the Ionian shelf (SE Sicily)[END_REF]. The observation of strongly thinned continental crust at the foot of the escarpment cannot be explained without invoking a major phase of rifting and crustal thinning in the Mesozoic. But the tectonic history of the ME is not the main subject of the Dellong et al., 2020 article and it is thus only briefly mentioned in the introduction and accordingly cites earlier work [START_REF] Argnani | Malta Escarpment fault zone offshore eastern Sicily: Plio-Quaternary tectonic evolution based on new multi-channel seismic data: Tectonics[END_REF][START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF]. In reply to comment key point one, it is unclear which article Argnani is referring to when he states (lines 27 -28): "that at least in one case previous literature was not adequately cited", since he does not specifically mention the article that was allegedly "not adequately cited".

Wide-angle seismic data alternate interpretation

Deep and asymmetrical sedimentary basin

Argnani questions the applicability of the rift basin interpretation, given a lack of expression in the surface morphology (line 76). He also comments that other authors have extended this interpretation to the entire area between the Alfeo Fault and Ionian Fault and to be possibly related to serpentine diapirs [START_REF] Polonia | A depositional model for seismo-turbidites in confined basins based on Ionian Sea deposits[END_REF]. To clarify this point, the "rifting" proposed in [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF] affects only the upper plate, the Peloritan -Calabrian basement block, which structurally represents the continental backstop of the Ionian-Tyrrhenian subduction. We never discussed nor implied it could affect the oceanic basement further south between the Alfeo and Ionian Faults, although indeed the Western lobe of the accretionary wedge is down-dropped more than 1000 m as a whole relative to the Eastern lobe and other publications propose this extension based on deep seismic data and geodetic models [START_REF] Polonia | A depositional model for seismo-turbidites in confined basins based on Ionian Sea deposits[END_REF]Mastrolembo Ventura et al., 2014). Indeed there is abundant evidence of NE-SW to E-W oriented extension in the straits of Messina area: from the pure normal faulting focal mechanism of the Messina 1908 M7.2 earthquake, to the extensional mechanisms of recent moderate magnitude earthquakes, to GPS data indicating up to 1mm/yr extension across the straits of Messina [START_REF] Palano | GPS velocity and strain fields in Sicily and southern Calabria, Italy: Updated geodetic constraints on tectonic block interaction in the central Mediterranean[END_REF]. And there is the observation of the 10-15 km deep basin of sediments above a thin (continental to transitional) crust discussed in earlier work [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF]. However, as Argnani correctly points out, the surface morphology is not characteristic of a large-scale continental rift basin. The explanation likely lies in the unique tectonic setting that this particular "rift" at depth, is overlain by 10-15 km of highly deformed, primarily folded, (by compressional to transpressional deformation) accretionary wedge sediments. In a sense as soon as the space (at depth) is made available, it is immediately filled by the overwhelmingly large quantity of actively deforming accretionary wedge sediments, present all around, and ready to fill the gap. This unusual sedimentary-tectonic history will not produce a classic syn-rift depositional sequence.

Alternate interpretation of the Calabrian Block

The study region exhibits a highly complex 3-D lithospheric structure and as more profiles became available, the original interpretation of the DY-P3 profile evolved. We agree that the lower portion of the crustal block, initially interpreted as Calabrian continental crust along profile DY-P3 in fact represents oceanic crust of the downgoing Ionian Sea slab (Figure 1; Comment key point 2). We note in the text: "While along DY-P3 no slab was modeled, along the profile DY-P4, the slab is clearly imaged by the data from the land stations [...] The Moho depth along model distance 80-120 km on profile DY-P3 (31 km) corresponds to the depth of the oceanic Moho along DY-P4; however, the backstop-slab interface was not detected along DY-P3 (Figure 9 and Figure S6 in Dellong et al., 2020;electronic supplements)". This new interpretation is due to the better data quality of the well-coupled landstations than the seabottom instruments. So we agree with some parts of the interpretation offered in this comment, which is different from [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF], but less so with Dellong et al., 2020.

Oceanic crustal thickness, Gravity models and Figure 2 of the comment

As to differences in crustal thickness between the oceanic crust imaged along DY-P4 and DY-P3, the top of the oceanic crust along the DY-P3 velocity model, was never shown in [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF] nor 2020 as it was not directly observed in the DY-P3 OBS data. If the author of the comment is referring to the oceanic crust location presented in the gravity models, it is worth mentioning that this specific "best fit" model was built to fit the freeair gravity anomaly and consequently does not reflect our latest interpretation of the area. The starting point of these models was the DY-P3 velocity model only, without the knowledge of the DY-P4 velocity model. The resulting interpreted oceanic slab depth coming from the gravity models was at around 30 km but with a high uncertainty. Later, the DY-P4 velocity model was produced and the oceanic slab depth was found around 25 km more precisely. Finally, in Figure 2 of the comment, the author is using an arbitrary iso-velocity line as "top of Ionian slab". This iso-velocity line does not correspond to any layers from the model, and does not correlate well with the DY-P4 velocity model (20 km depth instead of 25 km) which may explain the differences in crustal thickness observed by the author of the comment.

Differences with the interpretation of DY-P3 proposed by Argnani

The interpretation of the DY-P3 velocity model proposed by the author of the comment disagrees with the one proposed in the Dellong et al., 2020 on two major points: 1/ the interpretation of the green layer of intermediate velocities that was interpreted as the upper-crustal layer of the Calabrian block and is interpreted as a sedimentary layer by the author of the comment. And 2/ the Continent-Ocean Boundary (COB) location along the DY-P3 profile and the western edge of the oceanic slab along this same profile.

1/ In both articles [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF][START_REF] Maesano | Deformation and Fault Propagation at the Lateral Termination of a Subduction Zone: The Alfeo Fault System in the Calabrian Arc, Southern Italy[END_REF] this hypothesis was investigated, but a satisfactory answer with solid proof could not be achieved that this layer was an upper-crustal layer. We still prefer our final interpretation for at least two main reasons:

-Along the southern profile DY-P1, the intermediate velocity layer (4.5-4.8 km/s) is of sedimentary origin, as demonstrated by the CROP reflection seismic data section, along which a well-stratified layer is imaged in this depth. Unfortunately no such high-resolution MCS data exist for the northern profile DY-P3, however, seismic velocities are substantially higher than along DY-P1 (4.9-5.1 km/s) indicating a non-sedimentary origin of this layer here.

-Second, arrivals on the OBS from this layer differ between arrivals from the overlying layer, indicating more resemblance to the arrivals at the WSW of the profile. The lateral change takes place along a narrow region but is gradual and no abrupt contact was determined.

2/

It is worth noting that many of the regions where the author of the comment proposed a different interpretation from ours are in the deepest portions or at the ends of our velocity model, where ray coverage is sparse or absent and thus poorly constrained by the wide-angle seismic data. Thus, with the Dionysus velocity models presented in this study the location of the COB along the DY-P1 cannot be precisely determined. With this interpretation it makes sense that the thinned continental crust in the middle of the DY-P3 profile could belong to the ME as discussed in the [START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF] because of the overall structure of the ME is similar between DY-P1 and DY-P3.

"There may be no lithospheric tear fault in the Ionian Sea and south of Mount Etna".

Comment key point 3, Line 109 and the Conclusion lines 129-130

We strongly disagree with this, the crux of Argnani's comment. We stand by the interpretation given in the body of earlier work [START_REF] Gallais | Propagation of a lithospheric tear fault (STEP) through the western boundary of the Calabrian accretionary wedge offshore eastern Sicily (Southern Italy)[END_REF][START_REF] Gutscher | Tectonic expression of an active slab tear from high-resolution seismic and bathymetric data offshore Sicily (Ionian Sea)[END_REF]2017;[START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF] that from Siracusa southward the central and southern portions of the ME show no sign of modern activity. Argnani also agrees with this [START_REF] Argnani | Malta Escarpment fault zone offshore eastern Sicily: Plio-Quaternary tectonic evolution based on new multi-channel seismic data: Tectonics[END_REF]Argnani Comment, this issue). Also, according to Govers and Wortel, (2005) propagation of a STEP fault in the upper crust typically involves alongside vertical-axis structural rotations. However, the adjacent Hyblean Plateau is paleomagnetically unrotated [START_REF] Cifelli | An AMS, structural and paleomagnetic study of quaternary deformation in eastern Sicily[END_REF], which excludes the ME as a STEP fault. Despite other different interpretations existing on the exact location of the tear fault system [START_REF] Polonia | The Ionian and 205 Alfeo-Etna fault zones: New segments of an evolving plate boundary in the central Mediterranean Sea?[END_REF][START_REF] Barreca | Fault pattern and seismotectonic potential at the south-western edge of the Ionian Subduction system (southern Italy): New field and geophysical constraints[END_REF] we do all agree that the STEP is located in the Ionian Sea and not in the Tyrrhenian region. The AFS located 50 km eastwards (along profile DY-P1) matches all the expected characteristics of a crustal scale tear-fault, specifically a 50 km long elongated basin with a 500 -800 m thickness of syn-tectonic Neogene sediments, bounded by transtensional faults as imaged by highresolution seismic images and morpho-bathymetry of the seafloor [START_REF] Gutscher | Tectonic expression of an active slab tear from high-resolution seismic and bathymetric data offshore Sicily (Ionian Sea)[END_REF]2017). Moreover, these bounding faults show continuity with large-scale crustal heterogeneities imaged by deep seismic profiling [START_REF] Cernobori | Crustal image of the Ionian basin and its Calabrian margins[END_REF][START_REF] Nicolich | Crustal structure of the Ionian margin of Sicily: Etna volcano in the frame of regional evolution[END_REF][START_REF] Polonia | The Calabrian Arc subduction complex in the Ionian Sea: Regional architecture, active deformation, and seismic hazard[END_REF][START_REF] Gallais | Propagation of a lithospheric tear fault (STEP) through the western boundary of the Calabrian accretionary wedge offshore eastern Sicily (Southern Italy)[END_REF][START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF]. A recent seismic imaging study, based on a network of industry profiles, confirms the perfect correlation between the geomorphological expression of the transtensional structures (elongate basin, linear strike-slip to normal faults) and their deep expression as active faults affecting the basement [START_REF] Maesano | The Calabrian Arc: three-dimensional modeling of the subduction interface[END_REF]. Here a quote from their work regarding the identification of the Alfeo Fault system as the STEP fault: "Here we show the lateral continuity of the STEP fault system at depth for over 150 km length, and confirm its importance as a lithospheric structure" [START_REF] Maesano | The Calabrian Arc: three-dimensional modeling of the subduction interface[END_REF]. The latest work by this group [START_REF] Maesano | Deformation and Fault Propagation at the Lateral Termination of a Subduction Zone: The Alfeo Fault System in the Calabrian Arc, Southern Italy[END_REF] reinforces the same conclusions: "we confirmed the role of the AFS as a lithospheric tear", though they report on decoupling between shallow and deep fault segments, due to the effect of the accretionary wedge sediments. We have explained briefly here (in this paragraph) and in greater detail in earlier work [START_REF] Gallais | Propagation of a lithospheric tear fault (STEP) through the western boundary of the Calabrian accretionary wedge offshore eastern Sicily (Southern Italy)[END_REF][START_REF] Gutscher | Tectonic expression of an active slab tear from high-resolution seismic and bathymetric data offshore Sicily (Ionian Sea)[END_REF]2017;[START_REF] Dellong | Crustal Structure of the Ionian Basin and Eastern Sicily Margin: Results From a Wide-Angle Seismic Survey[END_REF] why the shallow and deeply-rooted tectonic activity of the Alfeo Fault System qualifies it as the best candidate for the present-day lithospheric tear fault. It is clearly distinct from, and located 30 -80 km east of the Malta Escarpment over most of its length (except in the Catania -Mt. Etna region, where the two structures intersect).

Figure 1 :

 1 Figure 1 : 3D view of the crossing point between DY-P3 and DY-P4 above the bathymetric map (Figure 1 & 2 of Dellong et al., 2020). The dotted line with question marks are showing a possible extension of the upper and lower crustal oceanic layers along the DY-P3 profile given the depth and location of the Oceanic crust along the DY-P4 profile. (Figure V.13 of Dellong PhD thesis, 2018).

Data availability statement

The ocean-bottom seismometer data used in this publication are accessible in standard Segy format upon request at http://doi.org/10.17882/52435