Chapter 1: Introduction

Forests management has always been a key issue and monitoring techniques are still evolving. Nowadays, it is important to be able to assess the amount, the diversity and the health of the forests on a territory. The main interest is economic for the raw resources they provide. The common goal is to estimate as accurately as possible their growth and yield. Furthermore, forests are an active part of our environment as they regulate the quality of water, soil and air. Therefore, new issues from the ecological perspective have risen lately when managing territories. The protection of the forests and of their biodiversity is a major concern for the biological research and the pharmaceutical industry as more diversity provides more basis elements for the development of new drugs. Forests management can also be an efficient tool to limit the current global warming as forests can store the atmospheric carbon dioxide and transform it into solid carbon. Furthermore, new bioengineering techniques suggest that it could be possible to modify the albedo of a forest canopy to increase the amount of reflected solar radiation in order to limit global warming.

The majority of monitoring techniques was developed at local scale for global estimations. From ground sampling to airborne data acquisition, the provided information related to local ecosystems always has to be extended to the whole territory. Scaling up the local data is used to produce estimations at a region-wide scale that cannot be directly assessed. Such field campaigns can require a lot of manpower, be time consuming. That is one of the reasons underpinning the development of the use of satellites for performing studies from space. Typically, remote sensing satellites use a variety of spectral bands to take pictures of the targeted forest (SPOT, PLEIADES, SENTINEL-2) as well as radar imagery (RADARSAT, ENVISAT, SENTINEL-1). These techniques are adapted to monitor changes and the health status of forests at a global scale, but they fail to assess accurately the density and height of forests.

LIDAR technology relies on the same principle as radar but uses electromagnetic waves of higher energy on the electromagnetic spectrum, typically from visible light to infrared bands. It is capable of travelling through thicker objects like tree crowns, thus giving an estimation of its density and depth. LIDAR has been mainly developed for ground (TLS) and airborne (ALS) applications. The use of LIDAR technology on satellites could then be a very powerful tool for forests management issues. Unfortunately, space environment is very challenging and developing new devices that can sustain the space conditions is very costly and time-consuming. LIDAR has only been used three times for space applications , by NASA. The first time, it was used on ICESAT from 2003 to 2009 in order to monitor the world ice caps. Recently, it was installed on ICESat2 and GEDI [START_REF]ICESat and GEDI missions[END_REF], respectively launched on September the 15 th and December the 5 th of 2018. ICESat2 will continue the missions of studying the ice cap and sea level evolutions started with the first ICESat. However, ICESat2 is using a new photon count technology for its sensors. Meanwhile, the main GEDI mission objective focuses on the monitoring of tropical and temperate forests and it has started to acquire data on March the 25 th of 2019.

The objectives of this internship take part in the LEAF-EXPEVAL project that aims to develop and validate simulation tools for the analysis of LIDAR signals on dense forests. As a supervisor of the project, the CNES would like to investigate the feasibility of a space LIDAR program for the study and monitoring of dense forests. In this context, the main objective of the LEAF-EXPEVAL project is to develop and validate tools for the simulation of a space LIDAR signal on forests [START_REF] Sylvie | LEAF -EXPEVAL proposition de recherché scientifique spatial pour[END_REF]. Once validated, the simulation approach will be used to give an insight on how the real space LIDAR should be calibrated before an actual launch could be planned. This project is conducted with the partnership of: . The proposed approach is to use the DART software for the simulation of LIDAR from space. The DART software is able to reproduce LIDAR signal emissions and receptions based on photons counting. The environmental conditions chosen for the simulations are dense tropical forests that are the most challenging environment for LIDAR technology. This is the reason why we chose this environment, because we want our LIDAR simulations to be the most possibly accurate under the most difficult conditions.

In this framework, the objective of my internship is to validate the approach developed for the simulation of large footprint LIDAR on tropical forest. To do so, I will compare real LIDAR signals acquired on a tropical forest plot with DART simulations realized using a detailed digital representation of the forest plot. The real LIDAR signals on tropical forest scene were obtained with an aerial LIDAR device (ALS), during the Paracou campaign of 2016 [3. The same scene was also scanned with a ground LIDAR device (TLS) in order to provide a 3D numerical mockup for the DART simulations.

The first section of this report focusses on Lidar signal and includes a brief presentation of Lidar technology and of DART, the software used for LIDAR simulations. In the second section I present the methodology developed to validate LIDAR simulations. Finally, I will present the results we obtained and the conclusions. The same principle as for radar is applied but LIDAR uses a monochromatic light source. A brief light pulse is emitted, reflects on a given target and is partly backscattered towards a sensor designed for the specific wavelength used. Usually, the laser pulse is Gaussian in time and space. When emitted, signal intensity is a 1D Gaussian with respect to time and the footprint intensity on the target describes a 2D Gaussian with respect to space. Those pulses are used to describe the given scene by recording the reflected signals, called waveforms. The interaction between light pulses and the target modifies the emitted laser signal. Depending on the optical properties (reflectance, transmittance) and on the position of the objects that make up the target, the return waveform intensity varies and can be used to describe the objects it has encountered. The light source and the sensor are usually located on the same device. As a consequence, there is a direct link between the time the waveform takes to return, the distance between the object it bounced off and the sensors. The relation between time and the distance to the target can be expressed as follow:

𝐷 = 𝑐 * ∆𝑡 2

Where D is the distance between the target and the sensor, ∆𝑡 is the time lapse between the emission and the reception of the signal and c is the speed of light. The position and orientation of the Lidar device are measured using a DGPS and an inertial unit. The position of the waveform echo (𝑥 , 𝑦 , 𝑧 ) can then be determined from the coordinates (𝑥 , 𝑦 , 𝑧 ) and the orientation (𝑋, 𝑌, 𝑍) of the sensor:

(𝑥 , 𝑦 , 𝑧 ) = (𝑥 , 𝑦 , 𝑧 ) + 𝐷 * (𝑋, 𝑌, 𝑍)

However, these equations are first order approximations, they only take account for lasers that reflect once on the target. When multipath occurs, it introduces noise onto the return waveform and delay the return signal. If those multiple reflections happen close to the ground or many times it can delay some signal echoes after the ground peak occurs. To solve this issue, it is commonly chosen to use a threshold on the intensity of the returning waveform. The part of photon reflections that comes from the ground forms a ground peak on the return waveform. This peak is weaken by vegetation density and in the case of very dense forests, like the tropical rain forest, it can even be completely absorbed by the vegetation. In this case, a threshold on the intensity can erase the ground peak. This is the reason why the tropical forest was used as a previous study case. Because this environment is the most challenging the LIDAR technology will encounter when mounted on the satellites. There are three types of LIDAR devices: terrestrial, aerial and spatial.

Terrestrial LIDAR, also called TLS, is mounted on a tripod that helps it stay level on any kind of terrain.

The device rotates to provide a 360° view of the surrounding, both vertically and horizontally. Due to the small divergence angle of the laser beam low the footprints, which are the areas lighten by each pulse, are small. Typically, the footprint width is at an order of centimeter at a hundred meters away from the device. This limits the possibility of multipath returns since the field of view for every pulse is very narrow. The maximum detection range depends on the transmitted power, target reflectivity, receiving aperture and detector sensitivity . Figure 2 shows what an indoor pattern of TLS pulses looks like, with the surrounding colored in blue and the LIDAR pulses colored in red: In the case of aerial LIDAR, called ALS, the carrier is usually a plane or a drone that flies above the area of interest. The height of the carrying device can go up to 2000 m. As a consequence, the footprints of ALS are larger than the ones of TLS. Typically, ALS footprints are ten to twenty centimeters wide. However, they remain very small for the study of a significant large area. Even for the study of a single tree, we need several measurements to fully represent it. Multipath may happens but the size of the sensors reduces its impact on the return waveform. Figure 3 shows how an ALS LIDAR is used to cover the target area below: The space LIDAR corresponds to the LIDAR devices mounted on satellites. On the low Earth orbit, the distance to the target can vary between 500 km and 1500 km. This means that the width of the space Lidar footprints is at the order of several meters (70 m for ICESat, [START_REF]ICESat and GEDI missions[END_REF]). The advantage is that it can cover very large area with a limited amount of measurements. The main drawback is that multipath, variation in ground level and atmospheric conditions introduce noises on the return waveforms. Since space LIDAR is a new technology, the optimal footprint shape and the required power level for the emitted pulses are still under investigation.

The LIDAR acquisition of the return waveform is a discrete signal. The intensity of the waveform is captured at a determined pace, called sampling time. The recording of the waveform can then be done under various formats. The two main formats used in this study are the Point Cloud (PC) and the Full WaveForm (FWF). With the FWF, the discrete samples of the waveform are stored and displayed as they were measured by the captor. Whereas for the PC, the discrete waveform signals undergo postprocessing to only retrieve the major intensity peaks. The different outputs of a LIDAR are presented in the following figure: Since a LIDAR tracks the light it itself emits, it is considered an active device of remote sensing instead of a passive one. It can operate on its own without any source of light from the surrounding. The advantage of such a technology is that it can be operated at any time of the day or night. The most commonly used wavelengths are located in the near infrared domain (NIR) between 700 and 1400 nanometers. In this spectral domain, the different vegetation types as well as the ground have a good reflectance coefficient. XXXX That is the reason why LIDAR opens up opportunities for the study of forest structure, including height and density measurements. However, the use of it in space is new and has not been completely proven to be effective yet. The latest GEDI mission provides few images such as the following on North American forests: This section has presented the LIDAR technology currently used. The next section will introduce the simulation framework for LIDAR using DART software.

II. 2. DART Description:

DART, Discrete Anisotropic Radiative Transfer, is a software tool for the simulation of radiative transfer in 3D. It has been developed by the CESBIO in Toulouse since 1993 and has been improved and supported ever since. It is currently used in many centers and Universities around the world such as NASA GSFC.

DART simulates radiation propagation for a wide range of electromagnetic waves, from visible light to thermal infrared (LWIR). Its radiation propagation modeling is iterative. It uses a flux tracking method ([10],), with a finite number of discrete directions for simulating the radiative budget and the images of optical airborne and satellite radiometers. For the simulation of terrestrial, airborne and satellite LIDAR signals (waveform and photon counting), it uses a so-called "RayCarlo" method that combines the Monte Carlo and ray tracking methods [START_REF] Gastellu-Etchegorry | Simulation of satellite, airborne and terrestrial LiDAR with DART (I): Waveform simulation with quasiMonte Carlo ray tracing[END_REF]. The environment into which the propagation takes place is called the scene and can be designed to represent a large variety of landscapes, such as trees and vegetation from natural forests or grassland and also buildings, concrete or glass elements from urban areas. The scene is modeled with cells, which amount depends on the location. Top of atmosphere cells consist of layers and the closer to the ground the location is, the more cells are used to describe the scene. The cells can then be filled with 3D elements, either turbid elements or a collection of planar elements. Planar elements are composed of triangles and parallelograms that are used to represent leaves and walls. They are defined by their orientation, size, and optical properties that can be chosen from a large collection such as Lambertian, Hapke, specular reflectance and isotropic or direct transmittance. These properties can also be implemented for a large collection of spectral bands. Simulated photons will interact with those planar elements upon meeting with them and depending on the optical properties. For example, elements describing objects with high transmittance, like glass, will let a lot of photons through. Turbid elements are a statistical representation of matter that either interacts with the photons or not. A turbid element has its properties defined on the cell it occupies. They are used for simulating fluids like air or water and vegetation, usually leaves and grass. When used for vegetation, the turbid elements are defined by an angular distribution and a volume density, representing infinitely small flat surfaces with optical properties such as Lambertian or specular. When used for fluids, the turbid elements are defined by the density of particles and by the particle properties: cross section, single scattering albedo and scattering phase function.

In this study, we focuss on the use of DART for the simulation of LIDAR data in aerial and spatial remote sensing conditions. The DART output of the LIDAR simulation can be under various formats. The ones we are interested in are the Point Cloud (PC) and the Full WaveForm (FWF) stored into LAS format [START_REF][END_REF].

The LAS format is a standard storage format used by the ASPRS (American Society for Photogrammetry and Remote Sensing). The same amount of bytes is used for each elements of the PC or the FWF. With the 1.4 version of LAS, the main information stored are: This section has presented the use of DART for LIDAR simulation. The two next chapters will present the data available for the LEAF-EXPEVAL project, the general workflow and the main steps already completed towards achieving the objectives of the project

 Xs,

Chapter 3 : Study site and Data

The ideal way to validate LIDAR simulations of large footprints would be to compare them with real spatial LIDAR signals on the same real forest target area. Unfortunately, such spatial LIDAR data is not available yet. This is why it was decided to validate the DART simulations for small footprint LIDAR data. The idea is to compare real aerial LIDAR data (ALS) and simulate ones on the same area of dense tropical forest. In order to run the DART simulation a 3D mockup representing the real forest scene needs to be built.. The scene also needs to respect the true position of the objects as well as the true terrain elevation, using a digital elevation model (DEM). The campaign to obtain the aerial and terrestrial LIDAR data sets was conducted in 2016 in Paracou, French Guyana

IV. 1. Calibration of the LIDAR device:

The LIDAR device used was a Riegl LMS-Q780 embarked on a plane. one important step was to calibrate the LIDAR device before the data could be compared with DART simulations. In order to do this, 14 different colored cloths and a spectralon were put on the ground to be used as targets for the calibration. The experiment consisted in different flight configurations over the different targets, with different reflectance. Both LIDAR and hyperspectral acquisitions were conducted at the same time. The following picture shows the targets used for calibration: The objective is then to solve the following LIDAR equation [START_REF] De Boissieu | LEAF-EXPEVAL: Calibration[END_REF]:

𝐼 (𝑡) = 𝐺 • 𝑆 • 𝑛 • 𝐷 • 𝑛 • 𝜌 • cos (𝛼) 4𝑅 • 𝑃
With:

 𝐼 : The intensity measured by the sensor (in numerical count)  𝐺 : The gain  𝑆 : sensibility of the sensor (A/W)  𝑅 : The distance to the target  𝑛 and 𝑛 : The respective transmittance of the system and the atmosphere  𝐷 : The opening diameter of the sensor  𝛼 : The incident angle of the ray on the target  𝑃 : The top peak power emitted during the pulse  𝜌 : The reflectance of the target Some of the parameters of the equation can be found in the manual of the LIDAR, such as 𝑃 , 𝑅, 𝛼 and 𝐷. The atmospheric transmittance is assumed to be equals to 1 for the flight passes at 450m. The different reflectance of the targets 𝜌 are known. The goal is then to identify the three remaining variables 𝐺, 𝑆 𝑎𝑛𝑑 𝑛 .

To simplify, these variables are grouped into the variable 𝐾:

𝐾 = 𝐺 • 𝑆 • 𝑛 • 𝐷
The LIDAR equation can be reshaped to fit 𝐾 in:

𝐾 • 𝜌 = 𝐼 • 4𝑅 𝑃 • 𝑛 • cos (𝛼)
From this last equation the parameter K can be identified using the different passes over the different targets. More details about the calibration step can be found in the calibration report from DeBoissieu [START_REF] De Boissieu | LEAF-EXPEVAL: Calibration[END_REF].

IV2 3D Mockup:

This section presents the construction of the mockup used for the simulation of LIDAR with DART. This mockup was built from the TLS data acquired during the Paracou campaign in French Guyana. The various steps for the construction this mockup can be described as follow :

1. Real TLS acquisition with geo-referencing, scans registration and filtering of the ground points 2. Filtering of wood and leaves echoes 3. Reconstruction of branches and trunks structures through a dedicated software 4. Reconstruction of the leaves density with a turbid medium (voxelisation) 5. Separation and Identification of trees 6. Analysis and implementation of the individual spectral properties 7. Building of the 3D mockup with both solid object structure for the trees branches and trunk and a turbid representation for the leaves density.

The first 1 to 5 steps of this list were done at the AMAP research unit. Most analysis T were done using the CloudCompare and Computree software. In the first step, the scene was scanned with an angular pace of 0.04°. In order to lower the point density, a 4mm minimum gap rule was applied.

CloudCompare was used to remove the repeating points within that 4mm gap. After this filtering the data were ready for the second step of classification wood/leave. The classification step was carried out using a comparative study of different separation mehods. The description of the different methods used can be found in [*]. The best method that was finally selected was a machine learning algorithm called Random-Forest. The classification into the two classes (wood and leave) was then validated and checked manually. In step 3, the tree trunks and branches are isolated from one another and their surface is mapped with 2D elements using the Computree software. The output of the step 4 is a .vox file. The voxels from AMAPVox can then be imported and linked with the cells onto the DART scene. After step 5, the output is a collection of objects in .obj and .ply format that we used to rebuilt the scene with the Blender software. The last two steps were done by the UMR TETIS team, see [START_REF] Vincent | Rapport CNES Leaf-ExpVal IRD AMAP[END_REF].

The attribution of the optical properties to the corresponding tree requires a 3D segmentation of the canopy vegetation. However, this operation is very complex in a dense forest and was not carried out. Instead, a 2D automatic segmentation was performed on the dominant trees of the canopy. To that aim the very high resolution multispectral image was segmented using the e-cognition software and the result was validated on the field.

Chap 5-Simulations and validation of the simulated signals

V.1 DART simulations

Indiquer juste ici que les simulations sont faites avec DART en paramétrant le lidar similaire au Riegl et en récupérant la trajectographie et les caractéristiques des pulses émis (tu peux éventuellement parler de Ptools4DARt mais risque de compliquer, à voir avec Florian Using the mockup that was previously constructed, we can launch the DART simulation of aerial LIDAR on our 3D mockup under the exact same conditions. For instance with the pulses we need to use the same wavelength, the same directions, the same intensity and the same duration.

Dire aussi sur quelles données le coefficient de calibration est appliqué pour rendre les données simu et réelles comparables

V.2 Development of a validation approach

The aim of this section is to provide a methodology for the validation of the simulation approach.

Once the small footprints are simulated, we need to aggregate both real and simulated footprint for comparison. Indeed, we cannot perfectly reconstruct the vegetation from the real forest scene. So, it is not possible to exactly reproduce every single waveform. We hope that the 3D reconstruction of the forest scene and the optical properties implemented are close enough to the reality, so that the aggregation of simulated and real waveforms will match to a larger scale. The aggregation enable us to estimate at which resolution the mockup we built is accurately representing the real forest scene. Our goal is then to compare the aggregations obtained from the real aerial LIDAR data and the DART simulations. We choose to perform the validation of this approach for both the simulations using full waveforms outputs (FWF) and the simulations using point clouds outputs (PC). In both cases, the real and simulated data sets will be compared on different criteria related to the aggregation of the aggregated waveforms or point clouds. The details of the aggregation methods and the computation of errors for comparison are developed in the following sections.

Aggregation methods:

The original waveform files were stored into .wdp format, so we extracted the waveforms and stored them into .las format file beforehand. Thanks to this operation, we can use similar algorithm for the treatment of FWF and PC data sets.

a. Full waveform aggregation method:

After the extraction, the full waveforms are stored in a dataframe. Each line of the dataframe corresponds to a single point on one of the waveforms and the column provides the following information:

 X, Y and Z :

The spatial coordinates of the point on the waveform  Intensity :

The intensity value of the signal at the waveform point  Waveform_ID : The identifier of the waveform that contains the point  Sample_ID :

The identifier of the point on the waveform  Sample_time : The time of acquisition between two sample points for this waveform We spatially subdivided the area of interest into "aggregation boxes". The size of these boxes depends on the resolution we want the aggregated signal to be. We defined the horizontal resolution of the boxes Lx and Ly and the vertical resolution Lz. The horizontal length of the boxes Lx and Ly is chosen to be equal and will vary from the actual size of the square area of interest (80m for Paracou) down to a length of 2m. The vertical length of the boxes will however remain constant at a length of 1m during the different case studies.

Once the domain is subdivided we can start to aggregate the waveform into the boxes. In reality, the return signal is a continuous intensity response but the sampling makes it discrete. The figure 10 shows what a typical waveform looks like while propagating through a subdivided domain: For every waveform, we want to add the contribution area to the corresponding box. The overall value in each box will then represent the aggregated intensity for that box. The discrete nature of the waveform signal means that we have to choose a method to approximate the contributions of the waveforms into the boxes. The method we choose for this approximation is to consider the contribution of each sampling point in its aggregation box completely, no matter how close to the box boundary it may be. Then we apply the rectangle rule. The intensity on a point is coupled with the sampling time of its acquisition to measure its integral. The figure 11 shows how the computation of the contribution for a waveform through a box is carried: The equation that calculated the aggregated intensity for each aggregation box is the following

𝐼𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 = 𝐼𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 (𝑥)
The algorithm that performed the computation of the waveforms aggregation can be described as follow:

For every box B: The methodology for the aggregation of point clouds is similar to the one for the FWF described above.


The scene is filled with point echoes that were aggregated into the same "aggregation boxes". This time we only need the spatial coordinates and the intensity of the echo. The method we used is simple: if the echo in located in box B, then we add its intensity contribution to that box B. This time we do not act on the values of intensity with sampling time because this was done during the construction of the echoes from the return waveform.

Comparison of simulated and real signal:

a-Data preparation

We compared the results from the DART simulation with the real airborne acquisition on the Paracou site in 2016. Both point clouds "PC" and full waveform "FWF" data set are investigated. The data sets are stored in .las format files. We rescale the vertical axis so that the height z=0 m always corresponds to the ground level. This rescaling may distortd the original waveform ray path depending on the steep of the ground below and its inclination from the vertical axis. However, since the aim of the study is to compare the "real" and the "simulated" data sets, it is not a big concern to apply the same rescaling to both data sets.

Here, we consider the comparison of "real" and "simulated" lidar data sets. We start by aggregating Lidar data into waveform (see section XXX) for different spatial resolution. The horizontal spatial resolutions Lx and Ly vary from the whole Paracou scene length, i.e. 80 meters, down to 2 meters. At first, the vertical resolution Lz is fixed at 1 meter height. The default setting conditions for the computation of the aggregated waveforms are sum up as follow:

-Horizontal resolutions take the following values: : Lx = Ly = { 80m, 40m, 20m, 10m, 5m, 2m} -Vertical resolution is Lz = 1m -There are no threshold on the vertical z axis (z is free) -There are no threshold on the intensity of the point as long as it is positive (I > 0)

In the following case studies, some of those parameters may vary alone or alongside others. For instance, the vertical resolution of the waveform aggregation will vary from 5m to 15cm (the length of a ns travel at lightspeed). The thresholds on intensity or height may vary as well. From the default setting parameters we also run a case study with a threshold of the Intensity of the echoes of the point cloud (I>7) as well as a case study with different vertical resolutions of the aggregation boxes (for ex. dz=2 m).

b-Comparison between real and simulated aggregated signals

To compare real and simulated signals, we choose to study different characteristics of the aggregated signal related to vegetation properties [START_REF][END_REF]::

-The intensity values of the aggregated waveforms or point clouds -The height of the peak of maximum vegetation intensity return -The heights of the Intensity quantiles from the ground level to the top of canopy (for all 5% quantiles)

The first two characteristics aim at evaluating the quality of the simulated signals with respect to intensity and the third set of metrics aims at evaluating its quality with respect to the shape. The comparison of these metrics between real and simulated signals therefore provides insight into the quality of the digital representation of a complex forest environment that was used as an input to the DART model.

To quantify the comparisons, we need to compute the approximation errors between real and simulated data set aggregations. For that matter, two types of error were considered: the classic L2 norm and the Normalized Root Mean Square Error (NRMSE). Those errors are computed to express the variation of the intensities aggregated at the box level, the height of vegetation peak and the height of the intensity quantiles.

The L2 norm of an error vector E is computed as follows:

‖𝐸‖ = ‖𝐼 𝐼 ‖ ‖𝐼 ‖ with ‖𝐼‖ = ∑ 𝐼[𝑥]²
And the NRMSE is defined by:

NRMSE = √𝑀𝑆𝐸 𝐼 with 𝑀𝑆𝐸 = ∑ (𝐼[𝑥] -𝐼[𝑥] )

c-Facteurs susceptibles d'induire des différences entre signaux réels et simulés Tu pourrais ici faire un point sur les facteurs pouvant donner lieu à des différences entre simulations et données réelles et indiquer comment tu vas les prendre en compte.

Ex : réflectance du sol non issue de mesures de terrain, représentation de la végétation sous forme de milieu turbide au lieu d'objet (ce point n'a pas été traité mais pourrait aussi expliquer le pic sol plus marqué), la sensibilité des détecteurs, et, pour le PC, la façon dont sont identifiés les échos) .

Pour la réflectance du sol et la sensibilité tu as proposé des moyens de détourner le pb (seuillage en h pour l'in, en intensité pour l'autre)

Chapter 6: Validation Results

In this chapter, we present the comparison of the real and simulated scene after aggregation for both the PC and the FWF data sets. The comparison and computation of errors is focused on the three metrics we identify previously: the aggregated intensity, the vegetation peak and the height of 5% intensity percentiles.

IV. 1. Full Waveform results:

In this section, we present the results obtained using the full waveforms data sets , coming from the real aerial LIDAR and the DART simulations. The results for default setting are presented for all investigated horizontal resolutions of the aggregation grid, from 80m (top left) down to 2m (bottom right) on figure 12: The real and simulated aggregated waveforms present similar shapes down to 10m resolution. As expected, the larger the resolution is, the better the real and simulated data fit. However, even when shapes are globally similar, the ground peak for simulated data is more pronounced than for real data. This means that ground spectral properties may have not been correctly set. Furthermore, the vegetation peak of maximum intensity seems to happen at the same height but with different intensity values. The simulated data set has a greater response, meaning that the vegetation optical properties may as well differ.

Figure 13 presents the aggregated waveforms at 80m (top left on figure 12) with the corresponding 5% height percentiles. At 80m resolution, the percentiles values fit correctly. However, we see that the shift in percentiles either comes from the initial greater ground peak for the simulated data or from the difference in intensity of the vegetation peak. To investigate the influence of the errors coming from the ground we used a height threshold. All points below 2m height were removed from the date sets. The figure 14 shows a collection of aggregated waveforms at different aggregation resolutions (80m/10m/2m) and with a threshold above 2m for the height of the measure waveform.

Figure14. real (red) and simulated (blue) aggregated waveforms with a ground threshold Without the ground peak, the real and simulated data sets seem to match better on the selected examples, but we need to confirm this with the computation of errors. The resolution limit for the mockup representation seems to be at around 10m as observed before.

The method of detection of the backscattering waveform may also be different between the real ALS LIDAR and the simulation. The algorithm and the intensity threshold from the constructor are unknown so we proposed to remove eventual noise by thresholding the intensity values below a given value.

The following figure show the results we obtained with a threshold on the intensity. The real data is much more affected when we apply an intensity threshold. The overall distribution of intensity on the real waveforms exhibit more low intensity than the simulated one. The overall difference at 80m resolution of energy (cumulated intensity) between the real and simulated data sets before thresholding was only of 2%. With the introduction of the intensity threshold the total difference of energy increases to 30%. There is a bias coming from the LIDAR device method for the acquisition of the real ALS. It seems that the detection threshold of the real LIDAR is lower than the simulated one..

The following figures 16 and 17 show different vertical resolutions for the aggregation of waveform. When we compare the red and purple curves to the blue one, we can see that the errors decrease proportionally with the vertical resolution as expected from the previous results. The thresholding of the ground point (green curve) leads to a small increase in similarity between real and simulated aggregated waveforms.,The optical properties of the ground could be improved, for instance, by lowering the ground reflectance. However, it seems that they can be neglected for the moment. The threshold on low intensity has worsen the comparison for all resolutions. This means that a lot of errors comes from the vegetation description. The optical properties of dense area vegetation need to be improved.

The following figures 19 and 20 show the errors computed on the others metrics: the vegetation peak and the height percentiles: When we look at estimation errors for vegetation peak and height percentiles of aggregated intensities, the DART simulations in full waveforms give good results down to 20m or 10m resolutions.

As expected, the main source of errors comes from the resolution and the limit of the voxelisation method seems to be at around 20m or 10m of horizontal resolution and 1m vertical resolution. From the intensity and ground thresholding, we see that some errors come from the representation of the ground and from the dense area of vegetation.

The following figures show the comparison of aggregated waveforms intensity by height layers with the default settings (figure 21) and with a threshold on the ground level (figure 22). In these figures, the green line represents the desired configuration where the plots from the simulated data are equal to the plots from the real one. Here again, the results from the comparison of vegetation peaks have a good match up to 10m resolution. On the smallest 5m and 2m resolutions we can see the impact of the ground threshold (see red circle on figure 23 and24). Figure 23 and 24 help us see that the greatest errors on the vegetation representation happen at low height. On the real ALS LIDAR, a lot of vegetation peaks were detected below 10m whereas they are not with the simulated data set (see bottom right on figure 24). That means that the lowest layers of the vegetation in Paracou were not well reconstruct on the 3d mockup.

On the other hand, in the case of the dominant trees at greater heights the comparison makes little mistakes even down to a 5m horizontal resolution.

Figures 25 and 26 show the quantiles comparison for the case study without any threshold (figure 25) and with a threshold on the ground level (figure 26): The same observations can be made with than with the vegetation peaks. The overall predictions of the percentiles are representative of the real scene down to 10 m. A lot of were located below 10 m for the real LIDAR whereas their corresponding height percentiles for the simulated data were located higher in the vegetation (see green zones). These errors at the lowest layers of vegetation may come from a bad representation of the vegetation structure and optical properties at these heights. We can also see the errors coming from the ground in the red zones on figure 25. These errors confirm that the ground properties should be adjusted, even though the overall errors is low. The NRMSE errors corresponding to those test cases are presented in the following figures 33 to 35. We did not edit the errors plots from the aggregation setting of 15cm vertical resolution because the error computed were too big to be compared with the other test cases. As a result, we conclude that it is not possible to correctly aggregated point clouds at a vertical resolution of 15cm. The Analysis of errors tells us that this time thresholding the ground points has better improved the comparison on both the percentiles precision and the vegetation peak position. This means that in this case, more errors come from the ground points proportionally as with the FWF data sets. With PC, the simulated vegetation structure is closer to the reality. Also this time, the threshold used to remove the low intensity points improves the comparison. This means that high intensity of backscattering zone are well represented. The fact that PC and FWF data sets exhibit opposite behaviors regarding the comparison of high intensity returns can makes us question the methods used during the acquisition of the data. In this case the mockup itself seems not to be the only source of errors. Still, the results show that the method gives a good approximation of reality down to 10m horizontal resolutions.

The following figures show the comparison between the settings with and without a threshold on the ground point. Figures 36 and37 show the aggregated intensity values at resolutions of 80m, 40m and 20m. From those more detailed comparisons, we can see that the PC simulations have a bias toward real intensity. However, from the point of view of aggregated intensity, the ground threshold have no impact.

Figures 38 and 39 show the percentile heights for all resolutions. In the case of the metrics on percentiles (figure 38 and 39) and vegetation peaks, the computed errors are lower when a threshold is applied to remove ground points, as expected from the computations of errors on figure 34. We can clearly see the disappearing of extreme points from percentiles comparison at the low resolutions of 10m, 5m and 2m (see red circles on figures 38 and 39). As it was observed before on FWF data sets, the errors at the lower level of vegetation are higher than the ones at greater height. However, this time, the errors at the low vegetation does not seem to be bias as they were with FWF: both simulation and real data detects vegetation in some areas where the other does not.

Chapter 5: Discussion of the results

Based on the previous observations we retain the following points:

 The simulation approach developed, i.e. running DART model on realistic tropical forest scenes obtained using TLS data and reflectance measurements to describe the structure and the spectral properties of the vegetation, respectively, has been validated down to 10m resolution:

The majority of errors came from the increase in the horizontal resolution. The aggregation of the small footprint waveforms at several resolutions has proved the mockup to be accurate and well representative down to 20m or 10m resolution, either with PC or FWF simulations. However, the representation of the vegetation is not adapted to obtain realistic waveform simulations for lidars with a footprint lower than 10 m.

 The vegetation structure representation can be improved: Both PC and FWF data, with the different analysis performed, show that the vegetation structure at low heights was producing errors. Both optical properties and vegetation structure could explain these errors coming from the low vegetation. Moreover, differences in the intensity of vegetation peaks at the canopy level show that the reflectance of the vegetation at high level could be improved as well. The contradicting results between PC and FWF may also indicate that the method to measure the echoes has introduced errors and could be corrected. New simulations should be carried out to adjust the turbid density used in the mockup. Also, a new approach using the repression with 3D triangles could be tested.

 The ground reflectance may be improved:

The errors coming from the ground peak returns can be neglected. Nevertheless, it shows that the optical properties of the ground could be improved, possibly by lower that reflectance in the DART simulation.

The perspectives for the LEAF-EXWPEVAL project are going to be divided on two axis:

o First, the approach proposed for the simulation of space LIDAR signals is valid for resolutions down to 20m at least. The DART simulation on the mockup can now be used to simulate space LIDAR signals and start the sensitivity studies for sizing a real LIDAR mission in space. o Second, additional investigations should be carried out to improve the simulations with DART.

The goal should be to improve the representation of the vegetation structure (especially for low altitude trees) as well as the ground reflectance. Different new analysis and data gathering could be considered: hyperspectral data from ALS for optical properties, adjusting the turbid density of the leaves and modifying the leaf angle distribution or using 3D triangles for the representation of the vegetation.

-

  IRSTEA (Institut National de Recherche en Sciences et Technologies pour l'Environnement et l'Agriculture) ; project leader. -CESBIO (Centre d'Études Spatiales de la Biosphère) -UMR AMAP (Unité Mixte de Recherche Botanique et Modélisation de l'Architecture des Plantes et de la végétation) -and CNES (Centre National d'Études Spatiales).

Chapter 2 :

 2 Background on Lidar signal II. 1. LIDAR technology: LIDAR (Light Detection and Ranging) is a remote sensing technology based on the echolocation of laser signals.

Figure 1 .

 1 Figure 1. LIDAR interaction principle (illustration adapted from [5])

Figure 2 .

 2 Figure 2. Principle of TLS (illustration from [6])

Figure 3 .

 3 Figure 3. Principle of ALS (illustration from [5])

Figure 4 .

 4 Figure 4. Two different types of LIDAR outputs (illustration adapted from [7])

Figure 5 .

 5 Figure 5. Vegetations reflectance: fully developed trees (green), small vegetation and dead leaves (brown) [8]

Figure 6 .

 6 Figure 6. GEDI first outputs, NASA[START_REF] Gedi | Image credits: Joshua Stevens (NASA Earth Observatory), Bryan Blair[END_REF]

Figure 7 .

 7 Figure 7. DART framework, illustration from DART user's manual[START_REF] User | CESBIO : Centre d'Etudes Spatiales de la BIOsphère -UMR 5126[END_REF]

  Ys, and Zs: The 3D position of sensor. Xs, Ys, and Zs values are used in conjunction with the scale values and the offset values to determine the coordinate for each point.  GPS Time: The GPS Time is the double floating point time tag value at which the point was acquired.  X(t), Y(t), Z(t): These parameters define a parametric line equation for extrapolating points along the associated waveform. The position (X,Y,Z) along the wave is given by: X = Xs + X(t) Y = Ys + Y(t) Z = Zs + Z(t)  Intensity: The intensity value is the integer representation of the pulse return magnitude.  Return Number: The Return Number is the pulse return number for a given output pulse. A given output laser pulse can have many returns, and they must be marked in sequence of return. The first return will have a Return Number of one.  Number of Returns: The Number of Returns is the total number of returns for a given pulse.  Byte offset to Waveform Packet Data: The waveform packets data are stored in the LAS file in an Extended Variable Length Record or in an auxiliary WPD file. The Byte Offset represents the location of the start of this LIDAR points' waveform packet within the waveform data variable length record (or external file) relative to the beginning of the Waveform Packet Data header.  Waveform packet size in bytes: The size, in bytes, of the waveform packet associated with this return. Each waveform can be of a different size.

III. 1 .

 1 Field Data . Décrire ici les données avec une partie terrain (forêt (incluant le TLS) et clairière avec bâches (avec mesures de réflectances) III.2. Airborne data Chapter 4 Workflow and preliminary stages Mettre ici le diagramme. Il faudrait rajouter la calibration radiométrique qque part sur ce diagramme ; c'est une étape avant la comparaison.Et dire que tu vas présenter la calibration et les maquettes utilisées en entrée des modèles, deux étapes importante pour faire les simulations et comparer données réelles et simulées ; Je mettrais d'ailleurs les maquettes avant la calibration qui intervient au moment de la comparaison.

Figure 8 .

 8 Figure 8.Calibration targets in Paracou 2016[START_REF] De Boissieu | LEAF-EXPEVAL: Calibration[END_REF]

Figure 9 .

 9 Figure 9. Validation Methodology illustration adapted from[START_REF] Mpili | Validation de la simulation du signal LiDAR sur la végétation forestière pour le dimensionnement d'un système LiDAR spatial[END_REF]

Figure 10 .

 10 Figure 10. Waveform through aggregation boxes

Figure 11 .

 11 Figure 11. FWF aggregation into one box

o

  Initialize a variable S for the aggregation of all waveform intensities in the box  Find every sampling point P that is located inside B  For every point P in B: Compute the variable 𝐼 = 𝐼𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 * ∆𝑇 o Add I to the sum S b. Point clouds aggregation method:

Figure 12 .

 12 Figure 12. Comparison of real (red) and simulated (blue) aggregated waveforms (from FWF) at different resolution (80m/40m/20m/10m/5m/2m)

Figure 13 .

 13 Figure 13. Aggregated waveform at 80m resolution with 5% percentiles (FWF)

Figure 15 .

 15 Figure 15 . real (red) and simulated (blue) aggregated waveforms with treshold on the intensity (I<7)

Figure 16 .

 16 Figure 16. real (red) and simulated (blue) aggregated waveforms at 15cm vertical resolution

Figure 18 .

 18 Figure 18. NRMSE errors of waveform aggregation at different horizontal resolution (from 80m to 2m)

Figure 19 .

 19 Figure 19. NRMSE errors of vegetation peak heights at different horizontal resolution (from 80m to 2m)

Figure 21 .

 21 Figure 21. Aggregated Intensity on the FWF data setsFigures 23 and 24 show the heights of simulated vegetation peak expressed with respect to real vegetation peak, with the default settings (figure23) and with a threshold on the ground level (figure24):

Figure 23 .

 23 Figure 23.Vegetation peaks on the FWF data sets

Figure 25 .

 25 Figure 25. Heights of 5% percentiles on the FWF data sets

Figure 29 .

 29 Figure 29. Real (red) and simulated (blue) aggregated intensity values with threshold on the intensity values (I<7)

Figure 33 .

 33 Figure 33. NRMSE errors on the aggregated intensity values (PC)

Figure 36 .

 36 Figure 36. Comparison of real and simulated aggregated intensity values (PC) at different resolution (80m/10m/2m)

Figure 38 .

 38 Figure 38. Heights of 5% percentiles on the PC data sets

IV. 2. Point cloud results:

In this section, we present the results obtained using the point clouds data sets . The data sets still come from the real aerial LIDAR and the DART simulations on the 3D mockup. The following figures 27 and 28 show the resulting waveform intensity from the aggregated point clouds, that we call aggregated intensity. Figure 27 shows the aggregated intensity at a resolution of 80m with the corresponding height values every 5% intensity percentiles: We observe on these PC data sets the same pattern as previously on the FWF data sets. The ground peaks are always present and more important on the simulated data than on the real ones. The waveform comparison looks visually good down to resolutions of 10m. However, this time, the vegetation returned intensities have greater values for the real data set.

The following figures 29 to 32 show aggregated intensity values for different settings:


José Pinho

Sabine Schmidt
email: sabine.schmidt@u-bordeaux.fr

Iris I Diallo

Hervé Derriennic

Hélène Fallou

Mario Lepage

Exploring the Susceptibility of Turbid Estuaries to Hypoxia as a Prerequisite to Designing a Pertinent Monitoring Strategy of Dissolved Oxygen

Keywords: estuary, Loire, Gironde, dissolved oxygen, hypoxia, turbidity, monitoring

Globally, there has been a decrease in dissolved oxygen in the oceans, that is more pronounced in coastal waters, resulting in more frequent hypoxia exposure for many marine animals. Managing hypoxia requires an understanding of the dynamics of dissolved oxygen (DO) where it occurs. The French coast facing the Bay of Biscay (N-E Atlantic Ocean) hosts at least a dozen tidal and turbid estuaries, but only the large estuaries of the Gironde and the Loire, are subject to a continuous monitoring. Here, we compared the DO dynamics in these two systems, in order to evaluate to what extent it is possible to predict the sensitivity to hypoxia of a tidal and turbid estuary based on its most common and easily accessible characteristics (the estuary surface area and length, liquid and solid discharge, suspended particulate load, water renewal time, presence of cities). Whereas the hyperturbid Gironde estuary is the longest (170 km), implying longer water and particle residence times, and is bordered by a large metropole (>750,000 inhabitants), only episodic summer hypoxia events (DO < 2 mg L -1 ) were recorded in its urban Garonne subestuary. In contrast, the turbid Loire estuary, smaller in terms of surface area and length, experiences permanent summer hypoxia. This demonstrates the inability to evaluate the susceptibility of a tidal and turbid estuary to hypoxia only by considering its general characteristics. We urge that there is a need for a pertinent DO survey strategy based on an initial assessment of DO in a specific estuary based on longitudinal investigations, in particular during the warm season (as higher water temperatures decrease DO solubility), for detecting the potential zones of minimum DO.

INTRODUCTION

The French Atlantic coast hosts numerous macrotidal estuaries that flow into the Bay of Biscay that are natural corridors for migratory fishes (Figure 1A). The two best known are those of the Gironde and the Loire, which are also the longest, at 180 km (considering the Gironde-Garonne system) and 100 km, respectively. However, there are also a dozen estuaries set geographically among them, of a smaller scale, such as those of the Charente (82 km), The range of the total river length (x) is given (the extension of the tidal section is not always defined for the small estuaries).

the Seudre (23 km) or the Auzance (5 km). The asymmetry of the tidal wave, when propagating upstream, creates a welldeveloped turbidity maximum zone (TMZ), with concentrations of suspended particulate matter (SPM) of several hundred milligrams of dry sediments per liter [START_REF] Toublanc | Formation and structure of the turbidity maximum in the macrotidal charente estuary (France): influence of fluvial and tidal forcing[END_REF]. This TMZ is defined as an estuarine region where SPM concentrations exceed 0.5 g L -1 in the Loire or 1 g L -1 in the Gironde [START_REF] Jalón-Rojas | Turbidity in the fluvial gironde estuary (S.-W. France) based on 10-year continuous monitoring: sensitivity to hydrological conditions[END_REF], 2017).

Extremely high SPM loads in transitional waters result in serious problems. In addition to silting, TMZ directly affects the water quality, in particular the oxygenation. First, the turbidity reduces the sources of dissolved oxygen (DO): it limits light penetration and thus photosynthetic activity [START_REF] Kromkamp | Possibility of net phytoplankton primary production in the turbid Schelde estuary (SW Netherlands)[END_REF], and gas exchange with the atmosphere [START_REF] Lanoux | Factors contributing to hypoxia in a highly turbid, macrotidal estuary (the Gironde, France)[END_REF]. In addition, the degradation of the organic content of the SPM consumes DO. This specific feature is likely to promote a chronic undersaturation of dissolved oxygen in turbid estuarine waters, with respect to the expected value in equilibrium with the atmosphere at a given water temperature and salinity [START_REF] Benson | The concentration and isotopic fractionation of oxygen dissolved in freshwater and seawater in equilibrium with the atmosphere[END_REF]. Severe depletion of dissolved oxygen down to borderline concentrations close to hypoxia (<2 mg L -1 ) is a major threat for most aquatic organisms [START_REF] Diaz | Spreading dead zones and consequences for marine ecosystems[END_REF][START_REF] Rabalais | Dynamics and distribution of natural and human-caused hypoxia[END_REF]. In estuarine waters, critical thresholds are even higher: here, we use a threshold of below 3 mg L -1 for hypoxia and between 3 and 5 mg L -1 for critical DO inducing stress for most organisms [START_REF] Curran | The oxygen requirements of a polluted estuary for the establishment of a migratory salmon, Salmo salar L., population[END_REF]. [START_REF] Delage | Oxythermal window drastically constraints the survival and development of European sturgeon early life phases[END_REF] showed that oxygen depletion induces sublethal effects at 70% O 2 saturation and lethal effects at 50% O 2 saturation during the early life phases of the European sturgeon Acipenser sturio. A minimum DO concentration of 5 mg L -1 is required all year long to preserve all resident and migratory fish species, including the most sensitive [START_REF] Breitburg | Effects of hypoxia, and the balance between hypoxia and enrichment, on coastal fishes and fisheries[END_REF]. Spatial extension of the deoxygenated zone is also critical for fishes, in particular for migrating species [START_REF] Alabaster | The dissolved oxygen and temperature requirements of Atlantic salmon, Salmo salar L., in the Thames Estuary[END_REF][START_REF] Woodland | Past and future habitat suitability for the hudson river population of shortnose sturgeon: a bioenergetic approach to modeling habitat suitability for an endangered species[END_REF][START_REF] Hrycik | Sub-lethal effects on fish provide insight into a biologically-relevant threshold of hypoxia[END_REF].

Dissolved oxygen is one of the key physico-chemical quality elements that must be considered under the EU Water Framework Directive (WFD) to reach a "good ecological status" of transitional waters [START_REF] Best | Dissolved oxygen as a physicalchemical supporting element in the water framework directive[END_REF][START_REF] Carvalho | Protecting and restoring Europe's waters: an analysis of the future development needs of the water framework directive[END_REF]. Hypoxia is indeed not the only stressor that impacts coastal ecosystems, but it also co-occurs with other factors (toxic substances, overfishing, urbanization) to decrease the ecological well-being of coastal waters (Committee on Environment and Natural Resources, 2010). Future global changes are expected to have negative effects on water quality: a warmer climate, by decreasing oxygen solubility and increasing respiration, may counteract the impacts of reduced organic pollution loading [START_REF] Bendtsen | Effects of global warming on hypoxia in the Baltic Sea-North Sea transition zone[END_REF][START_REF] Trolle | Effects of changes in land use and climate on aquatic ecosystems: coupling of models and decomposition of uncertainties[END_REF]. Based on the projected increases in temperature for 2050 and 2100, a 3D model of dissolved oxygen highlighted a trend toward a development of a permanent summer hypoxia in the urban Garonne subestuary [START_REF] Lajaunie-Salla | Future intensification of summer hypoxia in the tidal garonne river (SW France) simulated by a coupled hydro-sedimentary-biogeochemical model[END_REF]. There is then a need to record reliable observations in order to better understand DO dynamics and thus to improve preventive or remediation management strategies for hypoxia in transitional waters and to counterbalance the expected DO decrease under climate (temperature increase, salinization) and local (water abstraction, urbanization, waste water discharges) pressures [START_REF] Breitburg | Declining oxygen in the global ocean and coastal waters[END_REF]. However, there is a lack of systematic surveys in estuarine systems, and the small estuaries of a few km to 10 km are often neglected.

The objective of this work was to explore the possibility of assessing the susceptibility of tidal and turbid estuaries to hypoxia by considering their basic characteristics (estuary surface area and length, liquid and solid discharge, water renewal time). We have chosen as a case study the estuaries from the French Atlantic coast facing the Bay of Biscay. The strategy was to investigate the occurrence of hypoxia in the two best-known estuaries, the Gironde and Loire. The aim was to identify the potential prevailing conditions and to assess whether or not a survey of every estuary remains justified. First, we compared the main characteristics of the two systems likely to influence dissolved oxygen in order to forecast the oxygenation levels of their waters. Then, two different datasets (a long-term survey of the urban sector and longitudinal surveys) are presented in order to test the hypothesis. This study aimed to propose a strategy to test the hypoxic nature of an estuary as an indispensable prerequisite to further implementation of a DO survey.

MATERIALS AND METHODS

Study Site

The Gironde and Charente estuaries surround smaller-scale systems: from north to south the Boive, Falleron, Auzance, Vie, Curé, Lay, Sèvre-Niortaise, Charente, and Seudre (Figure 1). All of these estuaries that flow into the Bay undergo twice a day marine intrusion, with typical tidal ranges between 1.5 and 6 m at the mouth. Because of their macrotidal regime (tide height > 5 m) and the nature of their drainage basin, these systems present a maximum turbidity zone. The regional climate is oceanic: temperate with a warm summer. The perturbations, coming mainly from the Atlantic Ocean, bring rains almost all year round, but the rainfall peaks in winter whereas a summer minimum persists most often. Accordingly, fluvial discharges present a marked seasonal signal with a low-water period in summer and autumn. There is always the presence of, at least, a city, which is expected in 2019 to have wastewater collection and treatment systems according to the EU rules. However, urban wastewater treatment plants (WWTP) have variable efficiencies in treating organic matter (dissolved and particulate) and ammonia and in restraining storm sewage overflows.

High-Frequency Monitoring of Dissolved Oxygen

The Gironde and Loire estuaries have programs of long-term, high-frequency and multisite monitoring of essential physicochemical variables, called MAGEST and SYVEL, respectively. These networks measure the temperature, salinity, turbidity, and dissolved oxygen in the surface waters every 10 min. In 2004, the Gironde was the first to be equipped with a realtime measurement network; a detailed description is given in [START_REF] Etcheber | Monitoring water quality in estuarine environments: lessons from the MAGEST monitoring program in the Gironde fluvial-estuarine system[END_REF]. Four stations were constructed in the central and fluvial estuary, in particular in the urban Garonne subestuary in Bordeaux (Figure 1). Following this example, the Loire estuary was also equipped with automated stations at the end of 2006 [START_REF] Jalón-Rojas | Comparison of physical forcings affecting suspended sediments dynamics in two macrotidal, highlyturbid estuaries[END_REF]. The sensors were produced by Endress Hausser. The measuring range of the oxygen sensor (OXYMAX W COS41) was 0-20 mg L -1 (±0.2 mg L -1 ; [START_REF] Etcheber | Monitoring water quality in estuarine environments: lessons from the MAGEST monitoring program in the Gironde fluvial-estuarine system[END_REF]. The automated stations were found to be rather expensive to operate and were replaced in 2017-2018 by a new generation of autonomous multiparameter probes, called SAMBAT, developed by NKE Instrumentation (France). The measured parameters are temperature (-5 to +35 • C; ± <1%), salinity (range 0.1-42; ± <5%), turbidity (0-3000 NTU; ± <5%), dissolved oxygen (0-20 mg L -1 ; ± 0.1 mg L -1 ) and fluorescence (0-50 µg L -1 ; ± <5%). Data are transmitted daily to the dedicated data servers. This allows for a real time survey of water quality and of critical DO levels. It also permits the detection of potential problems (sensor drift) in addition to the normal required 3-month maintenance. The data are freely available on-demand 1,2 .

Longitudinal Surveys of DO Distribution

Fixed stations are very important for recording long-term observations with a high temporal resolution, but they cannot document the spatial extension of deoxygenation. Indeed, there are usually few instrumented locations and a low spatial resolution in a given system, as it is not possible to dramatically increase the number of stations due to the human and material costs. In this work, we applied a complementary strategy based on longitudinal surveys along the estuarine axis that we believe is the only way to determine the spatial extension of deoxygenation. In 2017 and 2018, field campaigns were realized in the Gironde and Loire estuaries. The transects were scheduled in order to begin the measurements at high tide from a site located upstream of an area where severe deoxygenation had already been reported. Then, we moved at a low speed downstream with a SAMBAT, maintained at 0.5 m below the surface, that collected a measurement every 2 min. The transects in the Gironde (Loire) estuary correspond to 172 (118) and 197 (90) measurements in 2017 and 2018, respectively. The boat was regularly stopped to send an optode (0-20 mg L -1 ; ±0.1 mg L -1 , SDOT, NKE Instrumentation) underwater in order to establish a vertical DO profile from the surface waters to the bottom. Simultaneously, water samples were collected to determine the SPM concentration (by filtering a known volume of water on a precombusted and preweighted GFF filter, error ± 2 mg L -1 ) in order to estimate the turbidity/SPM rating curve. In 2018, the particulate organic carbon (POC) fraction was also determined by direct combustion of these filters in an LECO CS 125 analyzer (error <3%). The data are available on an open-access repository [START_REF] Schmidt | Longitudinal Distribution of Salinity, Temperature, Turbidity and Dissolved Oxygen of Surface Waters in the Loire and Gironde Estuaries[END_REF].

RESULTS

Trend of DO in the Urban Sections

The estuaries of the Loire and the Gironde are characterized by the presence of a metropolis in their upstream reaches:

Nantes and Bordeaux, respectively (Figure 1). Hypoxia is often reported in the vicinities of such large cities, due to the input of organic matter and nutrients from wastewater discharges [START_REF] Garnier | Lower seine river and estuary (France) carbon and oxygen budget during low flow[END_REF][START_REF] Zhang | Factors contributing to hypoxia in the minjiang river estuary[END_REF]. The multiyear continuous record of dissolved oxygen in Nantes (Trentemoult site) and Bordeaux presents the same multiscale variability (Figure 2). First, the tide explains the short-term oscillations during the semidiurnal and fortnightly timescales.

The high currents at mid-tide, particularly during spring tides, enhance resuspension of particulate organic matter, which becomes available to heterotrophic bacteria and then promotes oxygen consumption [START_REF] Schmidt | Assessing and managing the risks of hypoxia in transitional waters: a case study in the tidal garonne river (South-West France)[END_REF]. However, the largest variability is clearly related to season. As the low salinity in both sites, typically below 2, is not expected to significantly influence oxygen solubility, the main controlling factors must be temperature, fluvial discharge or turbidity. The mean DO concentrations a 9.4 and 7.8 mg L -1 in the Nantes and Bordeaux surface waters, respectively. Half of the DO data in Nantes are rather high, from 8.0 to 10.9 mg L -1 , compared to that recorded in Bordeaux (5.8 to 9.7 mg L -1 ). The lowest DO concentration (0.74 mg L -1 ) was measured in Nantes (1.22 was the lowest in Bordeaux). In addition, it is noticeable that the Nantes surface waters experienced hypoxia events almost every year. Although the DO range is lower in Bordeaux, a single This trend is more obvious when considering the monthly mean (Figure 3). In Bordeaux, DO median concentrations are only slightly above the threshold of 5 mg L -1 in August and September. Nevertheless, at present, hypoxia events are still rare (Figure 3). In Nantes, the median DO is always above 6.7 mg L -1 , whatever the season. This might be seen as a conformable safety margin compared to Bordeaux. However, hypoxic DO concentrations were measured over a large temporal window, from May to November.

Longitudinal Surveys

The longitudinal surveys established in 2017 and 2018 correspond to highly contrasted situations. In 2017, it was conducted at the end of a heat wave and a dry first semester (mean fluvial discharge for the period from January to June: 596 and 493 m 3 s -1 for the Loire and the Gironde-Garonne, respectively). In 2018, the first half-year was particularly wet (1467 and 1173 m 3 s -1 ). This explains the large difference in salinity between the 2 years, especially for the Loire, as the transect ended almost at the mouth (Figure 4). The TMZ occurrence was also highly contrasted between the 2 years: during dry conditions (2017), the presence of the TMZ was very marked in the Gironde but to a lesser extent in the Loire. In 2018, the TMZ was nearly absent from the Loire and only recorded at the end of the transect in the Gironde. Fluorescence, used as a rough tracer of chlorophyll in turbid waters, was always higher in Loire waters and decreased with increasing SPM. Both the Gironde and Loire transects intercept the urban section. Nevertheless, it is clear that the minimal values along each transect are found elsewhere: systematically, approximately 20 km downstream from Nantes. For the Gironde-Garonne transect, it is much more variable: in 2017, the DO minimum was upstream from Bordeaux and associated with the highest temperature, and in 2018 it was downstream from the metropolis in the TMZ. The lowest surface DOs were always recorded in the Loire estuary. In addition, whereas the bottom and surface DO were rather comparable in the Gironde, the bottom Loire waters were clearly hypoxic in the deoxygenated zone.

DISCUSSION

Assumptions for Predicting Hypoxia Susceptibility in a Tidal Estuary

Several estuaries facing the Bay of Biscay are located between the Loire and the Gironde ones. Only these two large systems are instrumented for DO surveys and can provide the possibility of comparing theoretical and in situ assessments. The first objective of this work was to evaluate to what extent it is possible to predict the susceptibility of a tidal estuary to hypoxia. The theoretical evaluation relies on the comparison of the openaccess and common characteristics that are likely to influence oxygenation in a tidal estuary, besides turbidity (Table 1). The main forcing factors of DO dynamics in the Gironde estuary are temperature, fluvial discharge, particulate load and to a lesser extent tidal cycles and urban input (WWTP inputs, referring to metropolis inhabitants) [START_REF] Schmidt | Assessing and managing the risks of hypoxia in transitional waters: a case study in the tidal garonne river (South-West France)[END_REF]. We did not consider parameters such as excess nitrification, biological oxygen demand or organic matter and nutrient fluxes from the drainage basin, as this information is usually not available for small estuaries.

At the difference of Mediterranean-type estuaries, macrotidal systems, i.e., those with a tidal range >5 m, are well flushed and vertically partially or fully mixed, preventing a priori stratification [START_REF] Xing | A three-dimensional hydrodynamic and salinity transport model of estuarine circulation with an application to a macrotidal estuary[END_REF][START_REF] Warwick | Microtidal estuaries warrant special management measures that recognize their critical vulnerability to pollution and climate change[END_REF]. The surface and length of an estuary compared to the fluvial discharge determine the flushing efficiency and thus the water residence time, increases of which are adverse for dissolved oxygen. Supplies of organic matter and nutrients from the watershed and the cities are expected to increase with the size of the watershed surface area and the city population. Table 1 summarizes the main characteristics and their impact on DO, and highlights the most unfavorable case between the Loire and Gironde estuaries. The high residence times of water and of particles in the Gironde estuary are not favorable features for DO. Bordeaux is also a more populous metropolis compared to Nantes (Table 1; +23% in 2014). On the other hand, the Loire estuary is the receptacle of a larger watershed. The comparison of their different characteristics tends to describe the Gironde estuary as a site more exposed to hypoxia risk.

Hypoxia Occurrences in the Loire and Gironde Estuaries

The prediction based on their respective features was not confirmed by the DO dataset, whether this is assessed based on the long-term regional monitoring or the specific longitudinal surveys in 2017/2018. Over the past 15 years, the urban sections of the two considered estuaries have registered hypoxia events, though less frequently in the urban Garonne subestuary (Figure 3). In addition, the longitudinal surveys clearly demonstrate the presence of a permanent summer deoxygenation to levels below the hypoxia threshold downstream of Nantes. Along the Gironde estuary, the position of the low oxygen zone is more variable, and it can be up-or downstream of Bordeaux, depending on the fluvial discharge and the climatic conditions over the last few months. The localization of a large city (urban runoffs, WWTP discharges) does not appear to be one of the main criteria to locate the DO minimum along the estuarine axis.

Although the two estuaries are close spatially, there are differences in their DO dynamics. Table 2 summarizes the Pearson's correlation coefficients between daily averaged DO and various parameters (temperature, discharge, turbidity). There are notable differences for the same site for the correlation coefficients depending on the considered period (whole year or low-DO season) and among sites for a given period. On an annual basis, temperature is the main controlling factor of DO variability (-0.85) in Bordeaux waters, far beyond turbidity and discharge.

In the case of the Loire, at Nantes, turbidity is equivalent to temperature. The differences are even larger when considering only the periods during which low DOs have been measured (Figure 4); indeed, in the case of the Atlantic estuaries, the only common factor is the absence of a low DO during winter and early spring. From May to November, there is an obvious decoupling between the two estuaries. Whereas discharge and temperature have their main impact in the Gironde, SPM is the predominant controlling factor of deoxygenation in Nantes (-0.63). This was not expected considering the huge concentration, higher level up to 1 g L -1 , and extension, of several tens of km, of the TMZ in the Gironde, that limits exchanges with the atmosphere. One explanation could be the POC content of SPM. In the TMZ core, POC is on average 1.5%, with no seasonal variation, and 3.3% in the Gironde and Loire estuaries, respectively [START_REF] Etcheber | Particulate organic carbon in the estuarine turbidity maxima of the gironde, loire and seine estuaries: origin and lability[END_REF]. In 2018, the wet conditions permitted flushing the TMZ toward the mouth, and surface SPMs had a POC content ranging between 5 and 1.9% in the Gironde and 15 and 5.2% in the Loire (Figure 4). In the Loire, phytoplankton, a source of extremely labile organic matter, could contribute up to 50% of the total POC load due to eutrophication [START_REF] Etcheber | Particulate organic carbon in the estuarine turbidity maxima of the gironde, loire and seine estuaries: origin and lability[END_REF]. This is corroborated by the higher fluorescence content recorded during longitudinal surveys in the Loire. In the Gironde, the lowest POC of the 2018 longitudinal survey was associated with the TMZ presence (Figure 4), and in fact, the high turbidity of the waters prevents primary production during most of year and then limits the production of labile SPM, demonstrated by the low fluorescence concentration [START_REF] Irigoiena | Light limitation and distribution of chlorophyll pigments in a highly turbid estuary: the gironde (SW France)[END_REF]. This accounts for the differences in the two estuaries: the higher organic content of the particles enhances oxygen consumption in the Loire estuary. Following the categories defined by [START_REF] Kemp | Temporal responses of coastal hypoxia to nutrient loading and physical controls[END_REF], presentday hypoxia can be defined, depending on their occurrence; as persistent seasonal in the Loire estuary driven by a moderate TMZ during the warm season that does not impede primary production and subsequent organic matter degradation, and episodic in the Gironde estuary during peculiar hydrological and climatic conditions (i.e., droughts, heat waves).

CONCLUSION

Dissolved oxygen is one of the physicochemical quality elements that need to be considered to define the "good ecological status, according to the EU rules." The comparison of the Gironde and Loire estuaries demonstrates the difficulty in extrapolating the water quality of transitional waters only by considering basic and easily accessible characteristics. It is not possible to overcome preliminary acquisitions of DO concentrations of a tidal estuary to assess its susceptibility to hypoxia. The resulting question is then how to design a monitoring network, and in particular how to select the best sites to monitor. Our recommendations are to define first the local DO dynamics in relation to, at least, the discharge and SPM load over 2-3 years to investigate different hydrological contexts. In the case of the Atlantic estuaries facing the Bay of Biscay, the most critical period is the warm season, and the effort can then be concentrated on July to August. In the absence of historical DO data, it is necessary to conduct longitudinal surveys over the total estuarine length, and to establish evenly spaced depth profiles. The minimum parameters to determine are DO, temperature, salinity, SPM, and if possible, the fluorescence and organic content of particles. This investigation can help to define the zones of minimum DO. From there, it is possible to define the best sites to instrument for longterm surveys. When the MAGEST network was developed in 2004, the choice was made to monitor the city center, due to the presence of two WWTPs. The recent longitudinal surveys have revealed that there are in fact two intermittent low-DO zones, down-or upstream of Bordeaux, depending on the hydrology, justifying the future deployment of additional sites.

FIGURE 1 |

 1 FIGURE 1 | (A) Estuaries facing the Bay of Biscay between the Loire and the Gironde. (B) The Loire estuary. (C) The Gironde estuary. Gray areas indicate the main cities.The range of the total river length (x) is given (the extension of the tidal section is not always defined for the small estuaries).

FIGURE 2 |

 2 FIGURE 2 | Comparison of high-frequency DO records in the urban section of the Loire (A; 452,433 measurements) and Gironde (B; 506,167 measurements) estuaries. Gray areas underline DO concentrations corresponding to hypoxia stricto sensus.

FIGURE 4 |

 4 FIGURE 4 | Longitudinal surveys of dissolved oxygen and ancillary data (temperature, salinity, suspended particulate matter, fluorescence) in the Loire and Gironde estuaries in 2017 and 2018. POC content of SPM was determined only in 2018 (9 and 6 samples collected in the Gironde and Loire estuary, respectively). The distance from the mouth is normalized by the total length of the corresponding estuary. The measurement dates were the June 21, 2017 and July 5, 2018 in the Loire estuary; the June 27, 2017 and June 28, 2018 in the Gironde estuary.

TABLE 1 |

 1 Comparison of estuarine characteristics likely to influence positively ( ) or negatively ( ) water oxygenation when increasing.

	Impact		LOIRE GIRONDE
	on DO
	Tidal range	m	6.36	6.3
	Watershed surface area	km 2	117.5	81
	Estuary surface area	km 2	102	635
	Length	km	100	180
	Mean annual liquid discharge	m 3 s -1	860	680
	Water renewal time in summer	Day	20	>100
	Mean annual solid discharge	10 6 ton year -1	0.2-2	1.5-3
	Estuarine SPM stock	10 6 ton	1	5
	Metropolis inhabitants (2014)	nb	603,000	760,933
	Underlined values highlight the most unfavorable case between the Loire and the
	Gironde estuaries.

TABLE 2 |

 2 Pearson's correlation coefficients between daily average DO (mg L -1 ) and parameters, i.e., temperature, discharge and turbidity for the complete time-series (all seasons) and for the period from May to October.

	DO		Temperature	Discharge	SPM
	(mg L -1 )		( • C)	(m 3 s -1 )	(g L -1 )
	Gironde	Bordeaux
	All seasons	3578	-0.85	0.54	-0.38
	May-October	1916	-0.55	0.60	-0.25
	Loire	Nantes
	All seasons	3445	-0.67	0.47	-0.63
	May-October	1716	-0.20	0.35	-0.63

Frontiers in Marine Science | www.frontiersin.org

June 2019 | Volume 6 | Article 352

http://www.loire-estuaire.org/dif/do/network

http://rogir.epoc.u-bordeaux1.fr/Ayant_Droits/donnees/ Frontiers in Marine Science | www.frontiersin.org

ACKNOWLEDGMENTS

The authors thank the MAGEST and SYVEL network for the temperature and dissolved oxygen time series, and Dominique Poirier (EPOC) for the carbon determination. The authors sincerely acknowledged Mickaël Vallée, one of the last professional fishermen of the Loire estuary, for the local cruise completion.

FUNDING

This research was funded by the AFB, the French Office for Biodiversity (OFB) and by the OASU, (Aquitaine Observatory of the Sciences of the Universe). The authors warmly thank the long-term support of the Adour-Garonne water agency (AEAG) to the MAGEST monitoring network.

AUTHOR CONTRIBUTIONS

SS and ML conceived the research project. SS, ID, and ML collected field data during the longitudinal cruises. SS, ID, and HD operated the monitoring station at Bordeaux. SS analyzed the data and wrote the original draft. All co-authors commented on and approved the final manuscript.

Conflict of Interest Statement:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.