G Q Q Tran

T Dramais

J Van Emmerik

Némery

B Camenen
email: benoit.camenen@inrae.fr

N Gratiot

J.-A Cohard

F Gard

A.-T Nguyen

V Q Tran

G Dramais

T Van Emmerik

J Némery

Monitoring discharge in a tidal river using water level observations: application to the Saigon River, Vietnam

Keywords: Saigon River, water level, water discharge, tidal river, flood, stage-fall-discharge rating curve

Monitoring discharge in a tidal river using water level observations: Application to the Saigon River, Vietnam

Introduction

A good understanding of the hydrological cycle, and discharge in particular, in tidal rivers enables reliable forecasting and decision making by researchers and policy makers. Some priorities are generally put on the protections against floods, saline intrusion and the dynamics of pollutants because of the social, economic and political stakes they are linked to. Hydrological cycles in Low Elevation Coastal (and deltaic) Zone (LECZ, 0-10 masl) differ a lot from their upstream environments due to the tidal influence [START_REF] Van Driel | Vulnerability assessment of deltas in transboundary river basins[END_REF].

Although this is a strategic zone at the interface between land and ocean, physical and environmental data collected in this environment remain sparse because of inherent logistical difficulties due to the unsteadiness of the flow [START_REF] Taniguchi | Optimizing the waterenergy-food nexus in the Asia-Pacific Ring of Fire[END_REF]. On one hand, the river water discharge is influenced by tidal dynamics, which modulate hydrodynamics at high and low frequency and can reshape the geomorphology, with feedback loops on hydrodynamics [START_REF] Mao | Tides and tidal currents in the Pearl River Estuary[END_REF]. On the other hand, estuarine and deltaic zones are highly influenced by peak fresh water discharges, which are themselves modulated by tidal asymmetry [START_REF] Sassi | River flow controls on tides and tide-mean water level profiles in a tidal freshwater river[END_REF].

In the floodplain of the Saigon River, which is part of the LECZ of the south of Vietnam, vulnerability could be assessed according to four factors [START_REF] Mcgranahan | The rising tide: Assessing the risks of climate change and human settlements in low elevationcoastal zones[END_REF][START_REF] Van Driel | Vulnerability assessment of deltas in transboundary river basins[END_REF] : (i) Relative Sea Level Rise (RSLR) due to climate change and natural and anthropogenic subsidence;

(ii) Wetland ecological threat; (iii) Population pressure: and (iv) Delta governance (adaptivity, participation, fragmentation). The Saigon River flows through Ho Chi Minh City (HCMC), a highly populated area, where population density can reach up to 30,000 inhabitants/km 2 [START_REF] Nguyen | Nutrient dynamics and eutrophication assessment in the tropical riversystem of Saigon -Dongnai (southern Vietnam)[END_REF].

The unprecedented growth of HCMC induces pressure on the environment and especially water resources [START_REF] Van Leeuwen | The challenges of water governance in Ho Chi Minh City[END_REF].

A major issue in such regions is to measure the instantaneous discharge for a tidal river and to evaluate the residual discharge of the tide-affected river. Indeed, a simple stage-discharge relationship cannot be used because of the downstream influence of the tide. Stage-fall-discharge (SFD) rating curves are traditionally and successfully used to compute discharge at sites where variable backwater effects could affect a classical stage-discharge relationship [START_REF] Rantz | An empirical method of determining momentary discharge of tide-affected streams. Water-Supply Paper 1586[END_REF]. Such methods have recently improved thanks to Bayesian methods [START_REF] Petersen-Øverleir | Bayesian analysis of stagefall-discharge models for gauging stations affected by variable backwater[END_REF][START_REF] Mansanarez | Bayesian analysis of stage-fall-discharge rating curves and their uncertainties[END_REF]. However, the relation between stage, slope, and discharge has generally been considered too complex in tidal rivers to attempt to obtaining accurate discharges. [START_REF] Rantz | An empirical method of determining momentary discharge of tide-affected streams. Water-Supply Paper 1586[END_REF] proposed a graphical method using multiple correlation, which remains difficult to apply, especially when the tidal influence leads to discharge fluctuations as large as the residual discharge.

That is the reason why there exists nearly no hydrometric stations worldwide in such rivers. Recently, acoustic sensing techniques (through HADCP, Horizontal Acoustic Doppler Current Profiler) have been implemented successfully in some tidal rivers such as the Sacramento River, California [START_REF] Ruhl | Monitoring alternatives at the Sacramento River at Freeport, California: results ofthe 2002-2004 pilot study[END_REF] or in the rivers Berau and Mahakam in East Kalimantan, Indonesia [START_REF] Hoitink | Continuous measurements of discharge from a horizontal acoustic Doppler current profiler in a tidal river[END_REF][START_REF] Sassi | Discharge estimation from H-ADCP measurements in a tidal river subject to sidewall effects and a mobile bed[END_REF]. Such techniques generally correspond to a velocity index method. Similar method using an ADP (Acoustic Doppler Profiler) current meter was deployed on the Tanshui River, Taiwan [START_REF] Chen | Real-time discharge measurement in tidal streams by an index velocity[END_REF] or fixed ADCPs (Acoustic Doppler Current Profiler) on the Yangtze River at Xuliujing, China [START_REF] Zhao | Multiprofile discharge estimation in the tidal reach of Yangtze Estuary[END_REF][START_REF] Mei | Large addition of freshwater to the tidal reaches of the yangtze(changjiang) river[END_REF]. These methods remain however expensive (around 30keexcluding structural works and maintenance) and difficult to display in an urban place due to risk of vandalism.

In a recent work, [START_REF] Camenen | Estimation d'une courbe de tarage hauteur-dénivelée-débit pour une rivière influencée par la marée[END_REF] proposed to apply a stage-falldischarge rating curve to estimate instantaneous discharges in tidal rivers.

Such simple model can be calibrated using intense discharge campaigns achieved over a full tide cycle, which are easily available nowadays thanks to the ADCP technology. The first application of the method was based on field measurements made on the Saigon River in September 2016 and was particularly successful [START_REF] Camenen | Estimation d'une courbe de tarage hauteur-dénivelée-débit pour une rivière influencée par la marée[END_REF]. It presented some limitations for a very asymetric tidal wave (close to a tidal bore), which is not the case for the Saigon River.

In this paper, this simple low-cost water level-based discharge monitoring is applied the Saigon River for a two-year period. A calibration and validation of the method are first presented. The monitoring of two hydrological seasons from January 2017 to Dec 2019 is then discussed. Results obtained allow for examining and discussing two specific aspects : first, what is the hydrological pattern of the Saigon River at different time scales, from the event scale to years, and how is it influenced by the sea level ; second, how much the rainfall regime is directly influencing the Saigon River dynamics, in particular, how does the river react to extreme events? 2 Material and methods

Study site

The Saigon River is located in the South of Vietnam, in a low elevation coastal zone, i.e. between 0 and 10 m above mean sea level (MSL, Figure 1). The Saigon River takes its source in Cambodia and is 225 km long.

Its catchment area has a surface of 4717 km 2 [START_REF] Nguyen | Nutrient dynamics and eutrophication assessment in the tropical riversystem of Saigon -Dongnai (southern Vietnam)[END_REF]. The Saigon River is actually a complex river system, subject to several human and environmental interactions, including many canals while it crosses Ho Chi Minh City (HCMC) megalopolis, before flowing into the Dong Nai River and the coastal waters. Upstream the megalopolis, the river is regulated by the Dau Tieng dam, which was built during the 1980's, in order to mitigate saline intrusion and secure the fresh water supply uptake station of HCMC.

The situation of HCMC is all the more critical as 65% of the city is located at an altitude of 1.5 m above MSL [START_REF] Scussolini | Adaptation to sea level rise: Amultidisciplinary analysis for ho chiminh city, vietnam[END_REF][START_REF] Vachaud | Flood-related risks in ho chi minh city and ways of mitigation[END_REF].

The flow of the Saigon River is predominantly driven by tidal currents, which affects both the water level and water discharge, with regular exfiltration of water in some urban districts during high spring tides. Over the year, precipitation follow two contrasted seasons: the dry season, usually extending from November to April and the wet season from May to October, which gathers about 80% of the total rain [START_REF] Nguyen | Nutrient dynamics and eutrophication assessment in the tropical riversystem of Saigon -Dongnai (southern Vietnam)[END_REF]. Rain event can be severe and precipitation record can locally exceed 300 mm/day, even 400 mm/day. In this paper, we focus on one of the most severe rainy events that experienced the city. It occurred in November 2018, when Usagi typhoon hit HCMC. During this event, there was an increase of more than 300 mm in rainfall which caused flooding, material and human damages reported in several newspapers.

Water discharge estimation

The instantaneous water discharge was estimated by applying a stage-falldischarge (SFD) rating curve adapted from the general Manning-Strickler law, previously tested and validated by [START_REF] Camenen | Estimation d'une courbe de tarage hauteur-dénivelée-débit pour une rivière influencée par la marée[END_REF]. This model is based on the assumption of a pseudo-uniform flow and a prismatic section, which is valid for a low slope river with slowly varying water depth. Considering an almost constant flow between two water level measures, the Manning-Strickler equation can be applied : Another assumption made here is that the river section and reach are stable with no significant bed changes. It was verified by comparing the different bathymetries from both ADCP campaigns.

Q = KA w R 2/3 h √ S (
Two hydrological stations are needed to evaluate the hydraulic slope of the river: the water level z up and z dn were measured at Phu Cuong (upstream station) and Bach Dang then Thao Dien stations (downstream stations), which are named "PC", "BD", and "TD" hereafter (see Fig. 1a). Distance between stations is around L = 42 km (PC-BD) and L = 35 km (PC-TD), which is sufficient to have a significant difference in altitude and to estimate slope with a good resolution for this low slope system. Due to tides and flow oscillation, the slope oscillates between positive and negative values.

Also, the slope in Equation (1) (S = (z up -z dn )/L) does not correspond to the local slope at the point in which the discharge is evaluated, i.e. at Phu Cuong (PC). Indeed, because of the distance between the two stations, the flow discharge is slightly different at the two stations. However, this spatial offset can be translated as a lag time ∆t. As a consequence, the discharge estimation at PC can be written :

Q(t) = KA w (z up (t))R h (z up (t)) 2/3 |S(t + ∆t)| S(t + ∆t) |S(t + ∆t)| (2)
By specifying the river cross-section at PC (assuming it is relevant for the whole reach), one can easily evaluate A w and R h as a function of the water level z up . In Equation ( 2), two parameters need to be calibrated :

-K, the Strickler coefficient of the river reach supposed homogeneous between PC and BD, -∆t, a priori negative since the tide progresses from downstream to upstream.

The Strickler coefficient could vary depending on the discharge due to additional head losses at low flows (emergence of sills) or at high flows (flood-plain interaction). Nevertheless, such cases were not observed for the specific case of the Saigon River.

Water level measurement

To apply Equation (2), CTD Diver sensors were installed at PC and BD stations in September 2016, as reported in Figure 1. They measure pressure, conductivity and temperature every 10 minutes. Conductivity and temperature were considered as interesting proxy to interpret data series and were used to operate the post processing, but were not used directly in this study. stations and water levels measured by the SIHYMECC at the staff gauges for the selected period between 24th to 27th of June 2017.

Residual discharge

In order to study the net water discharge in the Saigon River and evaluate the hydrological cycle, the instantaneous water discharge needs to be averaged over four tide periods (T tide = 12.4 hours) to filter the 1st order semi-diurnal tidal signal from the data-series. As a consequence,

Q n (t 0 ) = t 0 +2T tide t 0 -2T tide Q(t)dt.
In the same way, the tidal averaged water level can be ob-

tained such as z ta (t 0 ) = t 0 +2T tide t 0 -2T tide z(t)dt.
The acquisition of data at high frequency (every 10 minutes) and over a longer period (two hydrological seasons) is a prerequisite to such evaluation.

3 Calibration and validation of the SFD rating curve

Calibration of the model using ADCP campaigns

Two Acoustic Doppler Current Profiler (ADCP) campaigns have been led

in September 2016 and March 2017 with a Rio Grande 600 kHz [START_REF] Dinehart | Repeated surveys by acoustic Doppler current profiler for flow and sediment dynamics in a tidal river[END_REF]. During 24 hours and for every hour, one gauging, i.e.

three transects, was realized at PC with a boat and a georeferenced ADCP mounted on it. ADCP campaigns were used to calibrate the water discharge estimation, calculated with Equation (2) [START_REF] Camenen | Estimation d'une courbe de tarage hauteur-dénivelée-débit pour une rivière influencée par la marée[END_REF]. K and ∆t were calibrated such as the modelled discharge fit to observed data (see also appendix A.1). Also, since it is very difficult to evaluate the exact vertical reference level of each station, these campaigns were used to optimize the measured difference z up -z dn using an additional parameter ∆z (i.e.

z up -z dn = (z up -z dn ) measured + ∆z)
To ensure the robustness of the hydraulic model, calibration were realized during both wet (September 2016) and dry (March 2017) seasons. These ADCP campaigns were also used to determine the river cross-section characteristics, to calculate A w and R h as a function of water level. Error in ADCP measurements were evaluated at 10% using a minimum error value of 100 m 3 /s since the conditions were adverse [START_REF] Le Coz | Estimating the uncertainty of streamgauging techniques using in situ collaborative interlaboratory experiments[END_REF]. We found the best results using a Manning-Strickler coefficient K = 26 m 1/3 /s and ∆t = -2.0 hours (see Appendix A.1). Results are presented in Fig. 3.

Very good results can be observed for the September 2016 campaign. Some slight error may be observed at the end of the first flow peak (Fig. 3a, at around 21 h, local time) but it is due to a pressure gauge outside water. For the 2017 campaign, results are not as accurate but still in good agreement with data. The clear asymmetric semi-diurnal tidal signal for this specific day may explain some of the differences [START_REF] Camenen | Estimation d'une courbe de tarage hauteur-dénivelée-débit pour une rivière influencée par la marée[END_REF].

Validation of the model using other discharge estimations

The monthly campaigns made by SIHYMECC also include some discharge estimation every hour for 48 h. They applied the velocity index method by measuring the water level and the depth-averaged velocity u index at one location y index [START_REF] Chen | Real-time discharge measurement in tidal streams by an index velocity[END_REF] : with α index = 0.8 a calibration coefficient, α index = U/u index where U = Q/A w the section-averaged velocity. Again, since conditions are adverse, we roughly estimated the error from this method equal to 15% plus a minimum error of 150 m 3 /s based on [START_REF] Ruhl | Computation of discharge using the index-streamflow: method in tidally affected areas[END_REF].

Q = α index u index A w (3) (a) (b)
Q [m 3 /s] Q (model) Q (ADCP)
The corresponding water discharge estimated from Equation ( 2) are in good agreement with data from SIHYMECC (Fig. 4). The model tends to yield smaller peak values for the inward discharge (negative), even if this trend is not observed for all tidal cycles. This tendency is more pronounced for the 2018 data (Fig. 4b). This could be the consequence of either a varying flow repartition throughout the river section during the ebb and flow not caught by the index velocity method or a bias on the slope estimation by the proposed model based on a slope downstream of the station.

(a) (b) In each panel, one can see that the tidal signal clearly predominates with a strong semi-diurnal harmonic. The cyclicity is particularly emphasized through the examination of a zoomed period in September 2018 (Figure 5b). The cyclicity is evidenced on both PC and TD stations, who present synchronous fluctuations, with a time laps of about 1 hour and 10 minutes, which corresponds to the duration of tide wave propagation between the two stations. The tidal magnitude is smaller at PC station (by 17%) than at TD station, which is physically consistent with quadratic friction dissipation: the further upstream the station is, the less it is influenced by tidal waves. Tidal range was measured to oscillate between -2.00 m and 1.50 m at both TD and PC stations.

A first glimpse of seasonality can be observed in the two-year time-series.

One can observe lower water levels in the rainy season (May to September) and higher values in the dry season (November to March). Semi-diurnal tidal forcing clearly prevails on the Saigon River. As illustrated in Figures 2 and5, the 14-days cycle presented both symmetric and asymmetric tides; a pattern which is particularly well developed in the estuaries located in the south of Vietnam [START_REF] Chen | Real-time discharge measurement in tidal streams by an index velocity[END_REF].

Continuous discharge estimation

By applying Equation ( 2), one can easily calculate the discharge at PC using water levels from both PC and TD stations. Water discharge is dominated by tidal semi-diurnal harmonics, as for water level data-series (Figure 6).

Instantaneous discharge ranged from -2000 to 2000 m 3 /s. The oscillation of the magnitude of tides, between spring and neap tides, is also well-identified (Figure 6b). Maxima appears at every syzygy between the sun, the moon and the earth. Overall, the analysis of instantaneous water discharge highlighted the predominant forcing of tides on any other forcing, despite contrasted and well-marked dry and rainy seasons.

Residual discharge

All data from the SIHYMECC field campaigns (which correspond to a 48 h experiment ≈ 4T tide ) were compared to the proposed model (Figure 7a). As expected, since the model yields smaller amplitudes (see Figure 4), one generally finds a good correlation but with few outliers, predominantly above the 1:1 line for negative discharge and bellow the 1:1 line for positive discharge (i.e. |Q model |lsim|Q(CEM )|). One should note also the many zero values for the SIHYMECC (CEM) data, which may be the consequence of a human bias (while reading the staff gauge) and partly explains uncertainties in experimental data and some of the scatter in the model response. A comparison of the net water discharge measured by the SIHYMECC and modelled by our group is reported in Figure 7b for the 18 SIHYMECC campaigns available.

One can realize the order of magnitude of difference between maximum instantaneous discharges and net discharges, which corresponds to the red box in Figure 7a. Despite the existence of few outliers, modelled points for net water discharge are well correlated with the estimation from measured data.

As expected, since the model yields smaller amplitudes, the net discharge estimations are generally smaller. The existence of outliers is not surprising looking at the sensitivity of such calculation while net discharge is one order of magnitude smaller than peak values (positive and negative). But surprisingly, while SIHYMECC estimations of the net discharge are always positive, the proposed model leads sometimes to negative values. Although difficult to explain, these negative values may result from a combination of factors, such a the asymmetrical tides and complex exchanges with groundwater and with the canal network. The regular occurrence of salt intrusions up to HCMC does confirm the possible occurrence of a net negative water discharge. Also, as discussed above, the model is very sensitive to a possible error in water level estimation a PC and BD, TD stations. Hence, for the net discharge, we evaluated the error equal to 20 m 3 /s per cm of error in the vertical axis (See Appendix A.3).

The net discharge time series over this two-year experiment is presented in The seasonality of precipitation in this tropical humid low lying area is expectedly a main driver of the annual hydrological cycle, but as shown previously, this driver is masked by the tidal forcing (Figure 8). The wet season is expectedly linked with a higher water discharge and a refill of the groundwater. At the opposite, the dry season is expectedly linked with a decrease of the net water discharge. This general pattern was broadly Monthly cumulated precipitation, tide-averaged water level along the Saigon and Dong Nai rivers and monthly averaged net-discharge are presented in Figure 9 for the 2017 and 2018 periods. 2017 and 2018 rainy seasons are well-marked, with an increase in precipitation from June to October. In 2018, precipitations are particularly high and extended in October and November. The tide-averaged water level follows the rain pattern for both increasing and decreasing phases, with a shift of two to three months.

Concerning the averaged water discharge, it was measured to fluctuate between -100 to 500 m 3 /s (no value is given for January 2018 since there was a too long gap in data), with a magnitude that was quite different for both hydrological years. In 2017, one can noticed that there exists an expected time lag between rain (first) and river discharge (few weeks later) but there is also some delayed and complex interactions with groundwater and/or water level.

The rainy season, which started in May, was followed by a rise in the net water discharge only one month later (June) and then, by an increase of the tide-averaged water level few weeks later (July, August). Also, one observed a drop in net discharge in August and September, while precipitation remains high followed by a rise in net discharge while precipitation dropped. It seems there are strong exchanges of water with groundwater and/or floodplain with a possible recharge from the groundwater and/or floodplain in May, August and September and a restitution to the river in November and December.

However, [START_REF] Van | Groundwater and river interaction parameter estimation in Saigon River, Vietnam[END_REF] Interestingly, these fluctuations are not correlated to precipitations and could partially explain the deficit in net water discharge compared to precipitation.

We suspect these fluctuations to be a consequence of the peak flow season of the Mekong River, which occurs in October in coastal floodplain and could influence sea water levels next to its river mouth (Chen et al., 2012a;[START_REF] Ho | Integrated urban flood risk management approach in contextof uncertainties: case study Ho Chi Minh city[END_REF][START_REF] Ha | Water and suspended sediment budgets in the Lower Mekong from high-frequency measurements (2009-2016[END_REF].

Effects of long-term precipitation

To get further in the analysis of the interaction between the Saigon River discharge and its surrounding floodplain, we examined a pluriannual monthly rainfall series (Figure 10). This series was associated with the severe Niño event in 2015-2016 [START_REF] Thirumalai | Extreme temperatures in Southeast Asia caused by El Niño and worsened by global warming[END_REF][START_REF] Ha | Water and suspended sediment budgets in the Lower Mekong from high-frequency measurements (2009-2016[END_REF]. In Figure The 2017-2018 and 2018-2019 hydrological seasons that are reported in Figure 9 support the hypothesis of a phase of recovery after these pluri-annual droughts. The 2017 wet season showed a priori an atypical hydrological response, with a normal rainy season from May to October, that is not in correlation with water discharge due to strong exchanges with groundwater.

There is indeed a shift between the discharge increase and the precipitation increase. This may result of years of deficit in the groundwater.

In 2018, the precipitation-river interaction was noticeably different. It seems the groundwater recovered from the drought leading to a more typical response of the river to precipitation. We may hypothesize that regional recovery of rain during the last wet seasons (2016)(2017)(2018) 2018), respectively. It clearly shows that floods observed in some districts of the city during the event were not directly driven by the overflow of water from the river, but certainly by intense run off over impervious streets. Also,

we should not underestimate the potential impact of the Dau Tieng dam management; large amount of water may have been retained during the flood and slowly reinjected in the river during the following weeks. Finally, since tide-average water levels in the Saigon River are relatively higher in late fall (Fig. 9), and because of the storm surge, evacuation of water may have been more difficult. This can be explained by the morphology of the delta, which is very flat and close to the sea shore. Results from [START_REF] Scussolini | Adaptation to sea level rise: Amultidisciplinary analysis for ho chiminh city, vietnam[END_REF] showed the influence of the sea level on HCMC flood risk. [START_REF] Ho | Integrated urban flood risk management approach in contextof uncertainties: case study Ho Chi Minh city[END_REF] argued that main hydrological explanations for the flood risk on HCMC area are the upstream floods including those of the Dong Nai and Mekong delta, local rainfall, land subsidence, and sea level rise. Our results indicate that local rainfall and sea water level may be the most impacting factors.

Conclusion

Tidal rivers and their floodplains are environments that are hardly monitored although they concentrate most of the human activities and form habitats for rich ecosystems. Flood risk is a particular issue for such environment and requires a good monitoring and a good understanding of its hydrological functioning.

Here, we applied a stage-fall-discharge rating curve adapted from the general Manning-Strickler law [START_REF] Camenen | Estimation d'une courbe de tarage hauteur-dénivelée-débit pour une rivière influencée par la marée[END_REF] to assess the instantaneous water discharge on the Saigon River. For the first time, the Saigon River dynamics was recorded at high frequency over two hydrological years, including the heaviest rain event ever recorded. Such tool appears to be a very efficient and low cost system for estimating discharge in tidal rivers as soon as the tidal wave is not too asymmetric and the river bed is stable (no general erosion nor aggradation). One issue remains the estimation of the slope since the model is sensitive to possible errors in water level measurements. This is particularly true for the Saigon River, which is tide dominated.

The analysis of the data-series highlights the driving role of tides on water mixing and water dynamics at hourly to monthly timescales. Once filtered from tidal effects, the averaged water level and water discharge points out a rather small contribution of run-off and a potential significant impact of flood plain and anthopogenic structures (Dau Tieng dam, canals, water pumping, etc.) to the hydrosystem response. The study of river response during Usagi typhoon illustrated how much the floodplain can attenuate the direct impact of intense precipitations and how critical is the sea water level to evacuate the flood. Also, the seasonal variability of the tidal-averaged water levels showed a possible influence of the nearby Mekong high flows that reduces the slope and limits the net discharge.

This study remains a preliminary study, which needs to be completed with some further studies on groundwater and flood plain dynamics, including the complex canal network in HCMC to get an accurate evaluation of the vulnerability of the megalopolis to flooding during extreme events and how this vulnerability will evolve during coming decades under a worrying context of subsidence and sea level rise.

A Sensitivity analysis

A.1 Error calculation on experimental data

A first sensitivity analysis was made to optimize main model parameters (Strickler coefficient K and the time shift ∆t) to experimental data. For both ADCP campaigns in September 2016 and March 2917 (see Figure 3), the estimation of the error was made using the mean square error :

E(Q) = 1 n Σ n i=1 [Q mod. (t i ) -Q meas. (t i )] 2 (4)
where n is the number of gauging data, Q mod. the modelled discharge and Q meas. the measured discharge, and t i the time of each of these gauging. As can be observed in Fig. 11, errors are significantly smaller for the first campaign for which the tide was more symmetrical. Indeed, E(Q) ≈ 200 m 3 /s for the September 2016 campaign whereas E(Q) ≈ 350 m 3 /s for the March 2017 campaign. From this sensitivity analysis, one would choose K = 25 m 1/3 /s and ∆t = -2.04 h. However, as can be seen in Figure 12a, using K = 25 m 1/3 /s would lead to some underestimation of the peak discharge values. Indeed, the error calculation is biased by the more numerous data points around the zero discharge value, which are slightly better predicted using a smaller value of the Strickler coefficient (see also Figure 3). Also, the effect of ∆t appears not as sensitive. It is also important to realize that water level data present a frequency of 10 min (1 min in September 2016), which does not allow to fit ∆t in such detail. Moreover, the position of the Thao Dien station (used for the 2017-2018 period) being a few km upstream of Bach Dang (used during the September 2016 experiment, see A.2 Sensitivity analysis on main parameters Equation (2) was calibrated by fitting two parameters: the Strickler coefficient K and the time shift ∆t. As can be observed in Figure 13, the Strickler coefficient has a linear effect on both instantaneous and net discharge. It is quite obvious since Q ∝ K. An error of approximately 8% in K would yield an error of 8% as well for both instantaneous and net discharge.

Since the coefficient ∆t leads to a time shift of the curves, it induces large errors around the flow reverse (Figure 14). However, the effect on the net discharge is very low; a shift of 12 min leads to an error of ±5 m 3 /s approximately.

A.3 Sensitivity analysis on input

The only input data of the model (Equation (2)) are the water levels measured at the downstream and upstream stations. The water level measured downstream will only affect the slope estimation whereas the water level upstream (at the reference station PC) will also affect the estimation of the wet section and hydraulic radius. In Figure 15 is tested an error of 1 cm in the downstream water level on the estimated discharge. The effect is inversely proportional to the discharge and errors can reach 300 m 3 /s for Q(t) ≈ 0 assuming an error of 1 cm. The asymmetrical shape of the figure is due the the time shift on the slope. On the other hand, error on the net discharge is quite limited but not negligible with an error of ±20 m 3 /s approximately.

In Figure 16 is presented the errors in discharge estimations for an error of 1 cm in the upstream water level. Surprisingly, effects are very similar (opposite) to the error in the downstream water level. Indeed, the sensitivity on the water slope is much higher than the sensitivity on the hydraulic radius or wet section.

Figure 1 :

 1 Figure 1: Location of the study site: the low elevation coastal zone of the Saigon River system, including the Saigon and Dong Nai rivers and Ho Chi Minh City (red points correspond to the position of pressure gauges).

  1) with Q the water discharge [m 3 /s], K the Manning-Strickler coefficient [m 1/3 /s], R h = A w /P w the hydraulic radius [m], A w the wet section [m 2 ], P w the wet perimeter [m], and S the energy slope [-] assumed equal to the water slope.

  PC sensor was immersed in the river bank of the Saigon River, close to Phu Cuong city, TD sensor was immersed in the heart of the megalopolis of Ho Chi Minh City, in Thao Dien district. Both sensors were downloaded at a bimonthly basis. Due to logistical constraints, significant risks of theft and of mechanical deterioration, TD station was displaced twice with a corresponding adjustment of parameters in Equation (2) (L and ∆t). The sensor was initially installed in Bach Dang, in the city centre in September 2016. It was moved to Boat House (BH), 8.5 km upstream of the Bach Dang point, from January 2017 to the 8th of March 2017. Then, since 15th March 2017, the sensor is located in Thao Dien Village (TD). Distance between those two last sites is only 0.9 km. To compensate water level measurements from atmospheric pressure fluctuations, a barometer was installed at the CARE center (Ho Chi Minh University of Technology, see Figure 1). Some drifts of the water level measurements were observed due to fine deposit accumulation leading to over-pressure. Water levels at PC and TD were thus corrected based on monthly campaigns made by Sub Institute of Hydrometeorology and Climate Change (SIHYMECC, also named CEM hereafter for Center of Environmental Monitoring) of Ho Chi Minh City. During these campaigns, water levels were measured every hour for three days at the staff gauge. In Figure 2, one can observed a very good agreement between CTD Diver data and SIHYMECC data (the water level H corresponds the the difference between the water surface level and a reference level of the station).

Figure 2 :

 2 Figure 2: Comparison between water levels H measured at PC (a) and TD (b)

Figure 3 :

 3 Figure 3: Water discharge Q comparison between model results and ADCP measurement campaigns in September 2016 (a) and March 2017 (b).

Figure 4 :

 4 Figure 4: Water discharge Q estimated at PC (blue line) in comparison to SI-HYMECC data for two contrasting hydrological periods in June 2017 (a) and in November 2018 (b).

Figure 5 :

 5 Figure 5: Two-year water level time series measured at PC and TD (a) with a zoom in September 2018 (b).

Figure 6 :

 6 Figure 6: Two-year water discharge time-series (a) with a zoom in September 2018 (b) at Phu Cuong (PC).

Figure 8 .

 8 Figure 8. Compared to monthly SIHYMECC campaigns, we observe more variability in the net discharge time-series. Some spikes are however not realistic and result from data gaps (January 2018, December 2019).

Figure 7 :

 7 Figure 7: Comparison of the instantaneous (a) and net (b) discharges obtained by the SIHYMECC (CEM) and the proposed model for all 2017 and 2018 SIHYMECC campaigns.

Figure 8 :

 8 Figure 8: Net discharge Q net estimate for the January 2017 -December 2019 period (full lines correspond to the model results while symbols correspond to values for the SIHYMECC campaigns using their own estimations -exp.-or using the model -mod.-; the grey line delimits a period of missing data at PC).

  estimated possible exchanges of water between the aquifer and the Saigon River of approximately 0.02 m 3 /s per km, which appears negligible. Effects of the large drainage and canal system around the Saigon River (i.e. of the Saigon floodplain) may prevail here. Also, anthropogenic influences such as Dau Tieng dam management or water supply pumping may significantly affect the net Saigon river water discharge. In 2018, cumulated rain also led to an increase of the net water discharge, but the rise of the Saigon River discharge was much faster and almost correlated to precipitation, which lets presume that soil saturation was higher and led to a direct response of the river to precipitation. Even if rain initiated late in the season, as compared with year 2017 (May and June in 2018 instead of April and May in 2017), the shape and cumulative rain during the rainy season was quite comparable to 2017. The net water discharge responded similarly to 2017 in June 2018. However, in 2018 the response to precipitation is stronger with a monthly averaged net discharge above 400 m 3 /s in July and August. Surprisingly, in October and November, while precipitation is particularly high in fall and superimposed with the extreme typhoon event in late November, one can observe a decrease of the net discharge. Significant floods have been observed during this period; a large part of the water volume from precipitation may have been spread over the flood plain.Looking at the tide-averaged water level (Figure9b) measured at Phu Cuong (PC), Thao Dien (TD) or Phu An (PA), but also at Vam Sat (station located downstream on the Dong Nai River at approximately 20 km from the see shore), a clear annual fluctuation can be observed with high water levels from October to April. Such high values yield an averaged negative slope leading some counter-effects on the net discharge (see Appendix A.3).

  10, one can observed a regular decrease of annual precipitation from 2009 to 2012, a minimum from 2013 to 2015 and an increase until 2019. The 2013-2015 period thus corresponded to a deficit of rain, which was evidenced by socio-political impacts described in newspapers, particularly during the dry season 2015-2016. This dry season was reported as the worst drought ever reported in 98 years of monitoring in the floodplains of Mekong and Saigon-Dong Nai hydrosystems.

Figure 9 :Figure 10 :

 910 Figure 9: Monthly rainfall P M (a), tide-averaged water levels H net (b) and monthly net discharge Q M,net (c) for the January 2017 -December 2019 period.

Figure 11 :

 11 Figure 11: Sensitivity of the Strickler coefficient K (a) and time shift ∆t (b) on the error calculation compared to ADCP campaigns.

Figure 12 :

 12 Figure 12: Comparison between the model and discharge measurement campaigns in September 2016 using K = 25 m 1/3 /s and ∆t = -2.00 h (a) and K = 26 m 1/3 /s and ∆t = -2.04 h(b).

Figure 13 :

 13 Figure 13: Sensitivity of the Strickler coefficient K on the instantaneous (a) and net (b) discharge.

Figure 14 :

 14 Figure 14: Sensitivity of the time shift ∆t on the instantaneous (a) and net (b) discharge.

Figure 15 :

 15 Figure 15: Sensitivity of the water level measured downstream on the instantaneous (a) and net (b) discharge.

Figure 16 :

 16 Figure 16: Sensitivity of the water level measured upstream on the instantaneous (a) and net (b) discharge.

Acknowledgements

This study was supported by the Rhône-Alpes region through the CMIRA Coopera financial support. We also wish to thank all the interns who worked on that topic throughout the last two years : Rémi Saillard, Lucas Barbieux, Roman Smoliakov and Swann Benaksas, as well as everyone involved in the two 24 hours field campaigns, which made the calibration possible. We are grateful to Dr. Nguyen Phuoc Dan from CARE center (HCMUT) and Dr Nguyen Van Hong from the sub-institute of hydrometeorology and climate change center (SIHYMECC) for providing us the field measurements of water level and water discharges.


Amandine Ramos

Jean-Patrice Robin

Lola Manizan

Cyril Audroin

Esther Rodriguez

Yvonne J M Kemp

Cédric Sueur
email: cédricsueur@iphc.cnrs.fr

Glucocorticoids of European bison in relation to their status: age, dominance, social centrality and leadership

Keywords: Bison bonasus, stress, sociality, collective decisions, ruminants, conservation

Stress is the body's response to cope with the environment and better survive.

While some studies suggest that it would be more stressful to be the dominant of the group, others support the opposite hypothesis. Several variables can actually affect this relationship, or even cancel it. This study therefore aims to make the link between social status and the basal level of stress of 14 wild European bison (Bison bonasus) living together. We collected faeces and measured the faecal glucocorticoid metabolites (FGM), a proxy of stress. We showed that FGM is linked to different variables of social status of European bison, specifically age, dominance rank, eigenvector centrality but also to several variable interactions. Hence, glucocorticoids seem to be a mediator of social life of bison. Preferential leaders in bison, i.e. the older and more dominant individuals which are also the more central ones, are less stressed compared to other group members. Measurement of such variables could thus be a valuable tool to follow and improve the conservation of species.

Introduction

Stress is the body's response to cope with the environment and better survive [START_REF] Boonstra | Reality as the leading cause of stress: rethinking the impact of chronic stress in nature[END_REF][START_REF] Crespi | Life history and the ecology of stress: how do glucocorticoid hormones influence life-history variation in animals?[END_REF]. For many social species, an individual's baseline stress level is generally associated with their dominance status [START_REF] Creel | Social dominance and stress hormones[END_REF]. However, while some studies suggest that it would be more stressful to be the dominant of the group [START_REF] Kotrschal | The relationship between social stress and dominance is seasonal in greylag geese[END_REF][START_REF] Mooring | Glucocorticoids of bison bulls in relation to social status[END_REF][START_REF] Sands | Social dominance, aggression and faecal glucocorticoid levels in a wild population of wolves, Canis lupus[END_REF], others support the opposite hypothesis [START_REF] Abbott | Are subordinates always stressed? a comparative analysis of rank differences in cortisol levels among primates[END_REF][START_REF] Mcleod | The relation between urinary cortisol levels and social behaviour in captive timber wolves[END_REF]. Several variables can actually affect this relationship, or even cancel it, such as the social system of the species, the seasonality of reproduction, domestication, personality, age or even the context of the study, i.e. in captivity or in the wild [START_REF] Creel | Social dominance and stress hormones[END_REF][START_REF] Engh | Patterns of alliance formation and postconflict aggression indicate spotted hyaenas recognize third-party relationships[END_REF][START_REF] Sapolsky | The influence of social hierarchy on primate health[END_REF][START_REF] York | The relationship of dominance, reproductive state and stress in female horses (Equus caballus)[END_REF].

The basal level of stress, measured via the secretion of glucocorticoid hormones, can also be intrinsically linked to the propensity to be a leader during group movements, i.e. to initiate the departure and be followed [START_REF] Ramos | How leadership could be used to manage domestic and wild ungulate herds[END_REF]. In equines, for example, the stallion usually moves behind the progression; this position allows it to maintain social cohesion by bringing its females together, i.e. herding behaviour [START_REF] Berger | Organizational systems and dominance in feral horses in the Grand Canyon[END_REF][START_REF] Feist | Behavior patterns and communication in feral horses[END_REF]. However, a study in Przewalski's horse species, Equus ferus przewalskii, showed that a change in the structure of the enclosure caused a significant increase in the secretion of glucocorticoids in the male. As a result, the latter then more often took the lead of collective movements [START_REF] Wolter | The role of an alpha animal in changing environmental conditions[END_REF].

A growing number of studies attempt to explain the phenomenon of leadership by focusing on individual characteristics (e.g. age, sex, physiological status, temperament, etc.) and on the social relationships that individuals maintain with each other (e.g. dominance, social centrality…). However, the majority only focus on the interaction between two or three factors (often age and dominance or dominance and sex) whilst these factors are often entangled and dependent on internal processes such as stress.

This study therefore aims to make the link between leadership, age, dominance, social centrality and the basal level of stress of individuals. Individual animals were identified and the stress level was measured by the concentration of faecal glucocorticoid metabolites (FGM, ng.g -1 ). We focused our study on semi free-roaming European bison (Bison bonasus), a wild ungulate living in groups with an average of 20 individuals [START_REF] Krasiński | Bisoniana XXIV. Free living European bisons[END_REF][START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF]. The only study reporting the interaction between glucocorticoid levels, individual characteristics and social structure in European bison was performed on captive herds and showed that cortisol levels were not influenced by sex or age (Klich et al., 2021a). However, the study site and the interaction between study site and visitor pressure were statistically significant. On the basis of behavioural observations made on leadership and sociality in this species [START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF](Ramos et al., , 2016a(Ramos et al., , 2016b[START_REF] Ramos | Collective decision making during group movements in European bison, Bison bonasus[END_REF] and the literature on the subject, we made the following hypothesis: the oldest females in the herd, which are the preferential leaders and the more dominant and central individuals, should have a higher average rate of FGM compared to the other group members.

Material and methods

Study area

Our study was conducted in the Kraansvlak area, a coastal dune area that it is part of the Zuid-Kennemerland National Park (3800 ha), the Netherlands (52° 23′ 16.1″ N, 4° 34′ 41.3″ E). This area was 330 ha large and consisted of six main habitats: sand dunes, grasslands, deciduous forest, pine forest, shrubbery areas and scattered water points. No natural predator of bison, i.e. grey wolf (Canis lupus) and brown bear (Ursus arctos), was present in the Kraansvlak area.

Kraansvlak is a nature core area within the National Parks where access to visitors is limited to guided excursions year round and a hiking trail accessible for the public during the months of September to March. The rest of the year the free-access walking trail remains closed. Further on only researchers and park rangers track the animals to monitor or study them.

Study Subjects

During the study period, the herd of European bison consisted of 14 individuals. Bison were identified at individual level (Table 1). The herd included eight adult females aged 4 to 15 years, two subadult females of 2 and 3 years old, two juvenile males both 10 months old and two juvenile females of 9 and 8 months old. The kin relationships between all individuals were known (Fig. 1).

Unlike numerous European bison populations living in natural reserves, the herd of Kraansvlak is not supplemented with food during winter or any other time during the year. Therefore, this herd shows the most natural foraging behaviour throughout the different habitats all year round. More details are given in [START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF] and [START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF].

Data Collection

Throughout the data acquisition period (from 5 April 2016 to 22 June 2016), we observed the bison herd from 09:30 to 16:30, for an average of 4 h per day depending on the climatic conditions as well as the time required to locate the group in the area. The total observation time for the entire study period is 156 h. We observed the bison at sufficient distance, minimum distance of 50 metres, which sometimes decreased slightly due to terrain geography or herd movement direction, to avoid setting off stress behaviours and ensure that our presence did not interfere with natural activity patterns and the welfare of the herd.

We recorded social relationships of individuals, dominance interactions, movement initiations (leadership) and collected faeces in order to obtain measure stress via glucocorticoids. Social and dominance data were collected using CyberTracker 3.0 Software [START_REF] Steventon | CyberTracker[END_REF], installed on a Samsung Galaxy Tab3 Lite. Leadership data was recorded with a HDR-CX320E Sony camera recorder.

Social relationships: Each 5 minutes, the identity of the nearest neighbour for each individual was recorded by the instantaneous sampling method [START_REF] Altmann | Observational Study of Behavior: Sampling Methods[END_REF]. When it was not possible for the observers to distinguish only one nearest neighbour (focal individual placed at the same distance of several individuals), the identities of these different neighbours were recorded. We used the eigenvector centrality as a measure of social centrality. This index was calculated using the nearestneighbour matrix implemented in SocProg 2.7 [START_REF] Whitehead | SOCPROG programs: analysing animal social structures[END_REF]. More details are given in Ramos et al. (2019[START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF].

Dominance interactions: All unidirectional agonistic interactions were recorded using the continuous sampling method [START_REF] Altmann | Observational Study of Behavior: Sampling Methods[END_REF]. For each interaction (see [START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF] for details), the identity of the emitter and the receiver were noted. Agonistic interactions were carried over an asymmetrical matrix to calculate dominance scores of individuals. This was analysed using SocProg 2.7 Software to measure the linearity of the hierarchy (h′ = 0.949, P < 0.0001, 10,000 permutations;

de [START_REF] De Vries | Measuring and testing the steepness of dominance hierarchies[END_REF] and to obtain a dominance index for each animal, the Modified David's Score (MDS). From these MDS, a dominance score, ranging from 1 (most dominant) to 14 (least dominant), was attributed to each bison. The dominance never changed during the observation period.

Movement initiations: We recorded the initiation phase of spontaneous collective movements with the continuous sampling method [START_REF] Whitehead | SOCPROG programs: analysing animal social structures[END_REF] by triggering the camera as soon as an individual moved away from the herd. The initiation of a movement was defined as the departure of the first individual (i.e. the initiator) out of the group diameter, for at least 20 steps in a constant direction (± 45°), with its head raised [START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF][START_REF] Ramos | Collective decision making during group movements in European bison, Bison bonasus[END_REF]. This head position has been described to be a good cue of individual motivation to change location [START_REF] Bourjade | Decision-Making in Przewalski Horses (Equus ferus przewalskii) is Driven by the Ecological Contexts of Collective Movements[END_REF][START_REF] Ramseyer | Decision making in group departures of sheep is a continuous process[END_REF]).

An initiation was finally considered to be a success when the initiator was followed by at least one congener [START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF][START_REF] Ramos | Collective decision making during group movements in European bison, Bison bonasus[END_REF]. If no bison followed the movement within the 15 min time window, the initiation was considered as a failure. The leadership score was then calculated by dividing the number of successful initiations per individual per the total of successful initiations (see [START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF] for the complete method).

Measures of faecal glucocorticoid metabolites (FGM)

The faeces samples were collected during the observations. For this, sterile 5 ml Eppendorf tubes were used, as well as disposable sampling spatulas. At least, one sample per day was taken, and all samples were identified on an individual level. The tubes thus filled were identified, dated and then immediately stored in a coolbag (containing an ice pack) in the field for a maximum period of 3 hours, then they were stored at -20 °C until laboratory analyses took place. At the time of sampling, we noted the environmental conditions (ambient temperature, sunshine or cloudiness, presence or absence of wind, rainfall level) and the time elapsed between defecation and sampling in order to determine whether fluctuations in environmental conditions and field constraints would have had a significant effect on the concentrations of FGM. To avoid contamination and deterioration, samples

were not taken on rainy days, when defecation took place in water or when individuals had diarrhoea. Indeed, after defecation, several factors, such as temperature, humidity, and other environmental conditions, may influence concentrations of FGM in the sample [START_REF] Millspaugh | Use of fecal glucocorticoid metabolite measures in conservation biology research: considerations for application and interpretation[END_REF]Möstl et al., 2002;[START_REF] Washburn | Effects of simulated environmental conditions on glucocorticoid metabolite measurements in white-tailed deer feces[END_REF]. Body condition might also confuse interpretation of FGM measures [START_REF] Heath | Body Condition and the Adrenal Stress Response in Captive American Kestrel Juveniles[END_REF]Smith et al., 1994). Because

bacterial enzymes are reported to increase or decrease levels of FGM if samples are not frozen shortly after sampling, only faeces dating from a maximum of 2 hours after defecation, were taken [START_REF] Millspaugh | Use of fecal glucocorticoid metabolite measures in conservation biology research: considerations for application and interpretation[END_REF]Möstl et al., 2002;[START_REF] Washburn | Effects of simulated environmental conditions on glucocorticoid metabolite measurements in white-tailed deer feces[END_REF]. Finally, only the faeces sufficiently isolated from each other and therefore reliably identified were collected. A total of 199 faecal samples were thus obtained and 103 were selected for analyses. The choice of samples to analyse was made in such a way as to have an equivalent quantity of samples per individual (N = 7±1), while making sure to have collection dates for each animal covering the entire period.

In collaboration with Professor Rupert Palme, the analysis of the level of FGM was carried out by the Physiology, Physiopathology and Experimental Endocrinology Unit of the University of Vienna (Austria). The protocol described in [START_REF] Palme | Steroid extraction: Get the best out of faecal samples[END_REF] and (Möstl et al., 2002) was used for the extraction and measurement of FGM using the EIA (Enzyme Immunoassays) technique. We obtained a value of FGM in ng/g-1 per sample (Table 1). This protocol was checked and approved covering many species including ungulates [START_REF] Christensen | Rider effects on horses' conflict behaviour, rein tension, physiological measures and rideability scores[END_REF][START_REF] Kovács | Salivary cortisol as a non-invasive approach to assess stress in dystocic dairy calves[END_REF][START_REF] Lemaître | Short-term telomere dynamics is associated with glucocorticoid levels in wild populations of roe deer[END_REF][START_REF] Touma | Measuring fecal glucocorticoid metabolites in mammals and birds: the importance of validation[END_REF][START_REF] Veissier | Heat stress in cows at pasture and benefit of shade in a temperate climate region[END_REF][START_REF] Wagner | Integration into the dairy cow herd: Long-term effects of mother contact during the first twelve weeks of life[END_REF]. We performed the measures twice, on wet and dry faeces and the values were correlated at 98.18 % (r = 0.9818, p < 0.0001). We chose to use FGM values on wet faeces as there is less protocol steps for measurement. The detailed protocol is given in the supplementary material.

Statistical Analyses

We used eigenvector centrality, leadership score and FGM per week during the observation period. Age and dominance rank were stable. Therefore, we could not implement individual identities as random factor in the next analyses as it was collinear with age and dominance rank (VIF > 4). Indeed, the potential collinearity between our predictor variables was tested by the calculation of variance inflation factor (VIF) from the R package car. In our different models (except Model 2, see below), age and dominance rank always showed a VIF > 4, meaning that they were highly correlated (adj r² = 0.79, F = 388, p < 0.0001, Figure 1a). Therefore, we decided to remove dominance rank as variable in the models 1, 3 and 4 as it is dependent on age, which is also shown in previous studies [START_REF] Bouissou | The social behaviour of cattle[END_REF][START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF][START_REF] Sosa | Impact of Group Management and Transfer on Individual Sociality in Highland Cattle (Bos taurus)[END_REF][START_REF] Sueur | Leadership linked to group composition in Highland cattle (Bos taurus): Implications for livestock management[END_REF]. Subsequently, multifactorial linear models (LM) were applied as follows:

- We tested leadership as influenced by the eigenvector centrality and not the converse, as usually done in studies on leadership [START_REF] Fischhoff | Social relationships and reproductive state influence leadership roles in movements of plains zebra, Equus burchellii[END_REF][START_REF] Jacobs | The influence of social organisation on leadership in brown lemurs (Eulemur fulvus fulvus) in a controlled environment[END_REF][START_REF] Van Belle | Collective group movement and leadership in wild black howler monkeys (Alouatta pigra)[END_REF], by using the function lm() from the R package car [START_REF] Fox | The car package[END_REF]. P-values for LM were calculated by Monte Carlo sampling with 10 000 permutations, using the function PermTest() of the R package pgirmess [START_REF] Giraudoux | pgirmess: Spatial analysis and data mining for field ecologists[END_REF]. Permutation tests for LM are well adjusted for moderate sample size and do not require normal distribution of model residuals [START_REF] Good | Permutation tests: a practical guide to resampling methods for testing hypotheses[END_REF]. We repeated three times the tests to check for the stability of the statistics (α = 0.05). All statistical analyses were performed with R, version 4.0.3 (R Development Core Team, 2009).

Results

Results of Model 1 (adj r² = 0.79, F = 67.29, p < 0.0001) show that the FGM is influenced by age (p < 0.0001, Figure 1b,c), eigenvector centrality (p = 0.002, Figure 1b), interaction Age*Eigenvector (p < 0.0001, Figure 1b) and interaction Age*Leadership (p = 0.0044, Figure 1c). Leadership and Leadership*Eigenvector did not influence the FGM (respectively p = 0.341 and p = 0.172).

Results of Model 2 (adj r² = 0.81, F = 150.2, p < 0.0001) indicate that the dominance rank is influenced by age of individuals (p < 0.0001, Figure 1a), by FGM (p = 0.006, Figure 1d) and interaction Age*FGM (p = 0.026).

Results of Model 3 (adj r² = 0.18, F = 8.75, p < 0.0001) show that the eigenvector centrality is influenced by FGM (p = 0.006), interaction Age*FGM (p < 0.0001) but not by age alone (p=0.408).

Finally, results of Model 4 (adj r² = 0.69, F = 45.94, p < 0.0001) indicate that the leadership score is influenced by age (p < 0.0001) but not by other factors or interactions (p > 0.353).

Discussion

In this study, we found that the FGM, a proxy of stress [START_REF] Chinnadurai | Validation of Fecal Glucocorticoid Metabolite Assays for South African Herbivores[END_REF]Möstl and Palme, 2002), is linked to different variables of social status of European bison, specifically age, dominance rank, eigenvector centrality but also to several variable interactions. Hence, glucocorticoids seem to be a mediator of social life of bison.

We had to remove dominance rank from three of our models, because this factor is highly determined by age in European bison [START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF] but also in other bovines [START_REF] Bouissou | The social behaviour of cattle[END_REF][START_REF] Sosa | Impact of Group Management and Transfer on Individual Sociality in Highland Cattle (Bos taurus)[END_REF][START_REF] Sueur | Leadership linked to group composition in Highland cattle (Bos taurus): Implications for livestock management[END_REF]. However, as shown in many studies, dominance rank is linked to glucocorticoid levels. While some studies suggest that it would be more stressful to be the dominant of the group [START_REF] Kotrschal | The relationship between social stress and dominance is seasonal in greylag geese[END_REF][START_REF] Mooring | Glucocorticoids of bison bulls in relation to social status[END_REF][START_REF] Sands | Social dominance, aggression and faecal glucocorticoid levels in a wild population of wolves, Canis lupus[END_REF], others support the opposite hypothesis [START_REF] Abbott | Are subordinates always stressed? a comparative analysis of rank differences in cortisol levels among primates[END_REF][START_REF] Mcleod | The relation between urinary cortisol levels and social behaviour in captive timber wolves[END_REF]. For instance, in captive animals, there are more aggression from the dominant towards the dominated for whom the inability to escape results in an increase of the glucocorticoid hormone secretions [START_REF] Creel | Social dominance and stress hormones[END_REF][START_REF] Kotrschal | The relationship between social stress and dominance is seasonal in greylag geese[END_REF]. Under natural conditions, a single study, carried out on olive baboons (Papio anubis), showed that the basal level of glucocorticoids was also higher in subordinate individuals [START_REF] Sapolsky | Endocrine aspects of social instability in the olive baboon (Papio anubis)[END_REF], while five other studies, grouping together several taxa, demonstrated a higher level of stress in dominant individuals [START_REF] Arnold | Reproductive suppression in male alpine marmots[END_REF][START_REF] Cavigelli | Behavioural patterns associated with faecal cortisol levels in free-ranging female ring-tailed lemurs, Lemur catta[END_REF][START_REF] Creel | Rank and reproduction in cooperatively breeding African wild dogs: behavioral and endocrine correlates[END_REF][START_REF] Creel | Social stress and dominance[END_REF][START_REF] Schoech | Reproductive Endocrinology and Mechanisms of Breeding Inhibition in Cooperatively Breeding Florida Scrub Jays (Aphelocoma C. Coerulescens)[END_REF]. Thus, according to the literature, even if there is no absolute consensus on the question [START_REF] Sapolsky | The influence of social hierarchy on primate health[END_REF], a high secretion of glucocorticoids is more generally observed in the dominant than in the subordinate individuals, in particular for species living in stable social groups [START_REF] Sands | Social dominance, aggression and faecal glucocorticoid levels in a wild population of wolves, Canis lupus[END_REF]. According to some authors, the high level of glucocorticoids in the dominant individual is not directly due to hierarchical status, but to agonistic behaviours expressed by the individual that allow him to access and maintain such a position [START_REF] De Villiers | Social dynamics and the cortisol response to immobilization stress of the African wild dog, Lycaon pictus[END_REF][START_REF] Sapolsky | Cortisol concentrations and the social significance of rank instability among wild baboons[END_REF]. In our study, we found that subordinate individuals have higher FGM levels than dominant ones but this is strongly mediated by age as subordinate individuals are younger individuals [START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF]. It is thus difficult to conclude about whether age or dominance is the prior factor affecting glucocorticoids, whilst studies have shown that, in same age groups of individuals, a correlation between dominance and FGM levels still existed [START_REF] Hasegawa | The effects of social exchange between two groups of lactating primiparous heifers on milk production, dominance order, behavior and adrenocortical response[END_REF][START_REF] Hubbard | Social dominance in beef cattle -A scoping review[END_REF]. Such group composition experiments would be interesting to test on bison, but this species being wild and vulnerable, it's not manageable like domestic groups.

In addition, we found a link between the FGM levels and the eigenvector centrality. More central individuals are less stressed. This social centrality is based on spatial proximity and thus is more a spatial centrality. Animals being more spatially central in the group can feel more protected [START_REF] Lopes | The influence of trees on the thermal environment and behaviour of grazing heifers in Brazilian Midwest[END_REF][START_REF] Mccann | Bed sites as thermal refuges for a coldadapted ungulate in summer[END_REF]. This correlation is also found with species showing grooming behaviour as central individuals are more groomed and this grooming decreases glucocorticoids level [START_REF] Brent | Social capital and physiological stress levels in free-ranging adult female rhesus macaques[END_REF]. Socially isolated individuals also show more stress (Sueur et al., 2021) as well in cattle [START_REF] Rushen | Opioid peptides and behavioral and physiological responses of dairy cows to social isolation in unfamiliar surroundings[END_REF][START_REF] Sec | Effect of different levels of social isolation on the acoustical characteristics of sheep vocalization[END_REF] as in other species [START_REF] Mumtaz | Neurobiology and consequences of social isolation stress in animal model-A comprehensive review[END_REF][START_REF] Sapolsky | Hypercortisolism Associated With Social Subordinance or Social Isolation Among Wild Baboons[END_REF] including humans [START_REF] Cacioppo | Social Relationships and Health: The Toxic Effects of Perceived Social Isolation[END_REF]. The relationship FGM*Eigenvector centrality seems to be mediated by the age of individuals. In fact, [START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF][START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF] showed that juveniles have similar centralities than adults because of their dependency but subadults have lower ones. The interaction here shows that FGM level decreases with centrality in juveniles and adults but not in subadults. The subadults are often meeting higher stress because they are less protected by their mother, especially if the mother has a new calf, need to find a place in the group and are subordinates.

A similar influence is found on the FGM level with the interaction Leadership*Age meaning that individuals who initiate more successful movements are less stress compared to other ones. This result can be explained by the greater experience and knowledge of older individuals which cope better with the uncertainty of changing location and resource availability [START_REF] Ramos | How leadership could be used to manage domestic and wild ungulate herds[END_REF]. The other explanation could be that individuals initiating movements do it because they feel the need to do so [START_REF] Fischhoff | Social relationships and reproductive state influence leadership roles in movements of plains zebra, Equus burchellii[END_REF][START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF][START_REF] Sueur | Predicting leadership using nutrient requirements and dominance rank of group members[END_REF] and by initiating these movements, they thus fill their needs and are less stressed. One may think that leaders are more stressed because of their position [START_REF] Harms | Leadership and stress: A meta-analytic review[END_REF][START_REF] Smith | Leadership and Stress[END_REF], but in humans, a study showed that leaders have lower levels of glucocorticoids probably due to a sense of control [START_REF] Sherman | Leadership is associated with lower levels of stress[END_REF].

Finally, and contrary to our initial hypothesis, we found that preferential leaders in bison, i.e.

the older and more dominant individuals which are also the more central ones [START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF][START_REF] Ramos | How can leadership processes in European bison be used to improve the management of free-roaming herds[END_REF], are less stressed compared to other group members. Other factors with a potential influence such as reproductive status, the transmission of social rank by previous generations, season (and consequently resource availability), or even parasitic infections should be studied further.

Conclusion

All these results are not only interesting for a better understanding of links between physiological and behavioural mechanisms in animal societies, but also for European bison management itself. Indeed, the species is listed on the IUCN Red List of Threatened Species and in need of receiving especially high attention in conservation and rewilding proceedings (Klich et al., 2021b;[START_REF] Pucek | European bison (Bison bonasus): Current state of the species and strategy for its conservation[END_REF][START_REF] Pucek | European bison: status survey and conservation action plan[END_REF][START_REF] Ramos | How leadership could be used to manage domestic and wild ungulate herds[END_REF]. Measurement of faecal cortisol metabolites, associated with other parameters, could thus be a valuable tool to follow and improve the adaptation and the welfare of herds in the semi-wild and in captivity (Klich et al., 2021a). This new knowledge could, for instance, allow us to take into account the importance of older females as group members as being stable and reassuring figures during captures, veterinary interventions or new reintroductions [START_REF] Ramos | How leadership could be used to manage domestic and wild ungulate herds[END_REF][START_REF] Ramos | The social network structure of a semi-free roaming European bison herd (Bison bonasus)[END_REF]. This could have a direct positive impact on survival of the herd.

  FGM~ Age*Leadership + Age*Eigenvector + Leadership*Eigenvector (Model 1) -Rank~ Age*FGM (Model 2) -Eigenvector~ Age*FGM (Model 3) -Leadership~ Age*FGM+ Eigenvector*FGM (Model 4)

Figure 1 :

 1 Figure 1: Main effects obtained between Faecal Glucocorticoid Metabolites (FGM, ng.g -1 ) and social status in European bison. a. Dominance in function of Age, b. FGM level in function of Eigenvector (arbitrary unit) and Age (months), c. FGM level in function of Leadership (relative frequency) score and Age, d. FGM level in function of Dominance rank.

Table 1 : Individual measures for the total observation period Individual Age (months) Eigenvector (A.U.) Dominance rank Leadership (relative ratio) FGM mean (ng.g-1)

 1

	Kareta (Kr)	180	0.29	2	0.23880597	272±89
	Kaga (Kg)	132	0.27	5	0.097014925	297±78
	Katarina (Kt)	132	0.32	3	0.253731343	210±52
	Katona (Ko)	102	0.32	1	0.02238806	187±47
	Ina (In)	78	0.17	6	0.104477612	177±70
	Moesja (Mo)	78	0.19	8	0.171641791	264±73
	Wisha (Wi)	78	0.22	4	0.052238806	227±47

Acknowledgements

We are very grateful to all the teams at Zuid-Kennemerland National Park, PWN, ARK Nature and all the partners of the Kraansvlak Bison Project for granting us access to the bison herd, for their valuable help during data collection and for their contributions during the various stages of the study. We thank Professor Rupert Palme and the Unit of Physiology, Pathophysiology, and Experimental Endocrinology from the University of Veterinary Medicine, Vienna for analysing the faeces cortisol data. Thanks also to Editage for reading and correcting the final version of this paper.