This thesis aims at better understanding the variability of the Red River plume in the Gulf of Tonkin (GOT) in the mid and far field area using a numerical modeling approach. Understanding the plume variability and the fate of the delta waters is of primary importance for an in-depth knowledge and a better prediction capacity of the ocean circulation and hydrology in the GOT, for an improved management of coastal waters and monitoring of the coastal ecosystems.

In this thesis, the SYMPHONIE numerical model is configured with a high-resolution variable grid and realistic forcings relying on the configuration of V. Piton (2019). It is run over 7 years, from 2010 to 2016 with daily outputs. Then, the model is assessed using several observational data and shows its good skill in most of the comparisons. The model configuration and verification provide the foundation for two subsequent studies below.

In the first study, the model results over a 6-year (2011-2016) period are used to examine the daily to interannual variability of the Red River plume and of three rivers whose mouths are nearby. The plume is identified with the help of the simulated passive tracers. The results show that the temporal evolution of the surface coverage of the plume is correlated with the runoff variations, but other processes, such as winds and tides, are involved. Using a K-means unsupervised machine learning algorithm, the main patterns of the plume and their evolution in time are classified in four clusters, analyzed and linked to different environmental conditions. In winter, the plume is narrow and sticks along the coast most of the time due to the downcoast current and northeasterly wind. In early summer, the southwest monsoon wind makes the plume flow offshore. The plume reaches its highest coverage in September after the peak of runoff. Vertically, the plume thickness also shows seasonal variations. In winter, the plume is mixed over the whole water depth, while in summer, the plume can be detached both from the bottom and the coast. The plume can deepen offshore in summer, due to strong wind (in May, June) or specifically due to a recurrent eddy occurring near 19°N (in August). Finally, the receiving basins of the different river systems in the GOT are identified, providing a general picture of the regions where the transport of materials from the river to the ocean is expected to have the greatest impact, for example in case of anthropogenic chemical substances leaked to the river. This first part was published in 2021.

The clustering analysis above shows that whatever the cluster, the plume is strongly affected by the wind. Therefore, the robustness of these results to the wind fields used to force the ocean model needs to be assessed. In the second study, I use an ensemble of simulations to assess the model response to perturbations added to the wind forcing. The sensitivity of the simulation presented in the first part is statistically evaluated by calculating the spread and the distribution of the variables of interest from an ensemble of 50 members. Due to computing and memory constraints, this study is performed over a short period, from June to August 2015, corresponding to the high runoff season. Firstly, the error of the forcing wind is estimated by comparing it with a satellite product. Then, its impact onto the model is assessed for surface and subsurface variables. For sea surface height, the average uncertainty is higher in the northern area, while for the sea surface temperature and salinity, it is higher near the Vietnamese coast and the Red River delta. Vertically, the uncertainty is highest at the surface for salinity and at the sub surface for temperature. The sensitivity of the river plume is then analyzed. The

spread of the plume area is highest in August, which is the same time as when the plume area reaches its peak. The clustering analysis shows some cluster shifts between different members of the ensemble, but the cluster that is most likely to occur is still the one from the reference simulation (with unperturbed wind). These limited changes suggest that the results of the cluster analysis of the reference simulation in the first study are indeed robust to the wind forcing errors. The uncertainty in plume thickness varies over time, typically less than 2m and can peak at 4m (for a total thickness of 10m).

Résumé

Cette thèse vise à mieux comprendre la variabilité du panache du Fleuve Rouge dans le Golfe du Tonkin (GOT) dans la zone proche de l'embouchure et plus au large, en utilisant la modélisation numérique. Comprendre la variabilité du panache et le devenir des eaux du delta est d'une importance capitale pour une connaissance approfondie et une meilleure capacité de prédiction de la circulation océanique et de l'hydrologie dans le GOT, ainsi que pour une meilleure gestion des eaux côtières et surveillance des écosystèmes côtiers.

Dans cette thèse, une configuration du modèle SYMPHONIE est mise en place avec des forçages réalistes et une grille à haute-résolution variable, sur la base de la configuration de V. [START_REF] Piton | Du Fleuve Rouge au Golfe du Tonkin: dynamique et transport sédimentaire le long du continuum estuaire-zone côtière[END_REF]. Il est exploité sur 7 ans, de 2010 à 2016 avec des sorties quotidiennes. Ensuite, le modèle est évalué à l'aide de plusieurs données d'observation et montre sa bonne habileté dans la plupart des comparaisons. La configuration et la vérification du modèle constituent la base de deux études ultérieures ci-dessous.

Dans la première étude, les résultats du simulation sur une période de 6 ans (2011)(2012)(2013)(2014)(2015)(2016) est réalisée pour étudier la variabilité journalière à interannuelle du panache du Fleuve Rouge et de trois rivières dont l'embouchure est voisine. Le panache est identifié à l'aide de traceurs passifs injectés dans la simulation. Les résultats montrent que l'évolution temporelle de la surface couverte par le panache est corrélée avec les variations du débit, mais que d'autres processus, tels que les vents et les marées, sont impliqués. En utilisant un algorithme d'apprentissage automatique non supervisé (K-means), les principaux régimes du panache et leur évolution dans le temps sont classifiés et analysés selon quatre clusters, puis liés à différentes conditions environnementales. En hiver, le panache est étroit et reste la plupart du temps le long de la côte en raison du courant côtier et du vent de nord-est. Au début de l'été, le vent de la mousson du sud-ouest fait s'écouler le panache vers le large. Le panache atteint sa plus grande couverture en septembre, après le pic du débit. Sur la verticale, l'épaisseur du panache montre également des variations saisonnières. En hiver, le panache est mélangé sur toute la colonne d'eau, alors qu'en été, le panache peut être détaché à la fois du fond et de la côte. Le panache peut s'approfondir au large en été, en raison de vents forts (en mai, juin) ou spécifiquement en raison d'un tourbillon récurrent se produisant près de 19°N (en août). Enfin, les bassins recevants des différents systèmes fluviaux dans le GOT sont identifiés, fournissant une image générale des régions où le transport de matériaux du fleuve vers l'océan sont susceptibles d'avoir le plus grand impact, par exemple dans le cas où des substances chimiques anthropogéniques seraient déversées dans le fleuve. Cette première partie a fait l'objet d'une publication en 2021.

L'analyse en clusters ci-dessus montre que, quel que soit le cluster, le panache est fortement affecté par le vent. Par conséquent, la robustesse de ces résultats aux champs de vent utilisés pour forcer le modèle océanique doit être évaluée. Dans la deuxième étude, j'utilise un ensemble de simulations pour évaluer la réponse du modèle aux perturbations ajoutées au vent forçant. La sensibilité de la simulation présentée dans la première partie est évaluée statistiquement en calculant la dispersion et la distribution des variables d'intérêt à partir d'un ensemble de 50 membres. En raison des contraintes de calcul et de mémoire, cette étude est réalisée sur une courte période, de juin à août 2015, correspondant à la saison de fort débit. Tout d'abord, l'erreur sur le vent forçant est estimée par comparaison avec un produit satellitaire.

Ensuite, son impact sur le modèle est évalué pour les variables de surface et de subsurface. Pour la hauteur de la surface de la mer, l'incertitude moyenne est la plus élevée dans le nord du bassin, tandis que pour la température et la salinité de surface, elle est plus élevée près de la côte vietnamienne et du delta du fleuve Rouge. Sur la verticale, l'incertitude est la plus forte à la surface pour la salinité et en sub-surface pour la température. J'analyse ensuite la sensibilité du panache de la rivière. La dispersion de la surface du panache est maximale en août, qui est aussi la période où la surface du panache est la plus grande. L'analyse en clusters montre quelques changements de clusters entre les différents membres de l'ensemble, mais le cluster le plus susceptible de se produire est toujours celui de la simulation de référence (avec le vent non perturbé). Ces changements limités suggèrent que les résultats de l'analyse clusters de la simulation de référence dans la première étude sont effectivement robustes aux erreurs de forçage du vent. L'incertitude sur l'épaisseur du panache varie dans le temps, elle est généralement inférieure à 2m et peut culminer à 4m (pour une épaisseur totale de 10m). Enfin, l'ensemble est vérifié en utilisant les jeux d'observations disponibles.

Acknowledgement

Doing a PhD is a special journey. My PhD is even more "special" because of COVID. The epidemic has affected us a lot from the beginning of 2020 (when I was a second-year student) until now (mid-2022). However, time has passed and we have learned that if we stand together, nothing can stop us from moving forward.

During this doctoral journey, I have received tremendous help and support from many people and organizations. First of all, I would like to express my deep gratitude to my supervisor, Pierre De Mey-Frémaux. This PhD could not be started without your help. You guided me from the beginning, you supported me when I applied for the fellowship, and you shaped my study.

I want to express my particular and deepest gratitude to my co-supervisor, Nadia Ayoub, who supported me day and night. She was the person to follow my every step along my PhD journey, not only at work but also in my social life. I will remember all the memories I had with you and your family. It was a pleasure for me to be your student.

I am extremely grateful to Thanh-Ngo Duc, my co-supervisor. Discussions with you about the atmospheric conditions helped me a lot to strengthen the method used in my study. Thank you for all your valuable advice not only during a year in Vietnam but also when I was in France.

I am also indebted to the other members of the research team. A big thank you to Patrick Marsaleix for your help in setting up the models and for all the problems related to the models. Also, you are the one who suggested me to try the passive tracer in the study. I am also grateful to Florence Toublanc for all the discussions throughout the thesis. Finally, thank you for introducing me to Python 5 years ago. I could not have undertaken this journey without Sylvain Ouillon, who initiated the link between me and my supervisor. I believe that the connections you have made have proven to be strong and undeniable! I would like to express my sincere thanks to the people on the thesis committee, Marine Herrmann and Anna Rubio. Although we only met once a year, I appreciated your advice and suggestions. I think your recommendations shaped the method I used in this study.

A big thank you to my friends in the lab, Violaine, Simon, Julia, Adélaide, Arne, Quentin, Micael, ... Thanks for the time we spent together. I will miss the international meal very much! Words cannot express my gratitude to my Vietnamese friends in Toulouse: anh Việt, anh Đa, Bích, Hạnh, chị Huệ, anh Thái, chị Ngọc, Thái con. The party we had almost every week has colored my social life. I will remember the memories we have together. I am also grateful to the people of USTH: chị Trà, anh Sơn, anh Bình, anh Hiền, anh Thường, anh Tuấn, anh Long, anh Quân, anh Kiên admin. I spent an unforgettable year with you in Vietnam. Special thanks to my friends, who I could always count on when I was stressed: Tuấn, Đức Anh, Hải dớ, Thu vb, Hoàng híp.

The last word, I want to give it to my family. My mother and father didn't have the chance to go to college. They started working before they were 18 and they have continued to work until now. They worked tirelessly so that I could have enough money to go to college. No words can express my gratitude for their sacrifices.

I would also like to thank the French Embassy in Vietnam and IRD (ARTS program) for funding this thesis. This work used the high performance computing resources of CALMIP (University of Toulouse III Paul Sabatier) and HILO (Hanoi University of Science and Technology). Thanks to the authors (Manh Cuong Tran, Florence Lyard, Sylvain Ouillon) who provided the observational data to validate the model.

List of figures

List of acronyms

I.1. This work

The study of river plumes is of great interest to the coastal oceanographic community because of the significant impact that plumes can have not only on local hydrology near the mouth, but also on ocean dynamics and ecosystems across the shelf. They also carry natural and anthropogenic contaminants that enter the ocean along with river waters ; the monitoring of river plumes is thus essential for water quality management policies. All of these issues are both more critical and intensified by the fact that coastal areas near river mouths (through estuaries, deltas, lagoons) are generally highly populated: anthropogenic pressure is higher on inland and coastal waters, while water quality and ecosystem health issues impact a large population. This is particularly true in tropical and subtropical deltas where sea level rise, land subsidence, flooding, and salinity intrusion are among the major threats in the context of global change (e.g. [START_REF] Szabo | Population dynamics, delta vulnerability and environmental change: comparison of the Mekong, Ganges-Brahmaputra and Amazon delta regions[END_REF].

The Gulf of Tonkin is bathed by the waters of the Red River, which reach the ocean after passing through a complex network of channels and distributaries forming the Red River delta. The apex of the delta is located about 50km upstream from the capital Ha Noi. The delta is home to a population density of about 1000 inhabitants/km 2 , and represents an area of intense agricultural and industrial activities. Understanding the physical and biogeochemical functioning of the complex system constituted by the Red River delta and its plume is essential to better face the challenges posed by global change.

In this context, the main objective of this work is to provide a better understanding of the variability of the Red River (RR) plume in the Gulf of Tonkin (GOT) at different scales.

This thesis work takes place within the framework of the LOTUS project. LOTUS is a joint Vietnam-France project, which aims to better understand the coastal ocean dynamics under the influence of the atmospheric variability and of the delta. This is a crucial knowledge to evaluate natural hazards impacts, to address environmental issues and to set-up sustainable coastal water monitoring systems. The work done within LOTUS involves Vietnamese researchers from several universities and institutions (USTH, HUS, IMER, IO, …) as well as French researchers from IRD, CNRS and French universities. In the framework of the LOTUS project, several studies in the GOT have been conducted. In addition to my work, two other PhD studies deal with the ocean dynamics in the GOT. In the PhD thesis of V. [START_REF] Piton | Du Fleuve Rouge au Golfe du Tonkin: dynamique et transport sédimentaire le long du continuum estuaire-zone côtière[END_REF] at LEGOS (Toulouse, France), a high-resolution configuration in the Gulf of Tonkin and downstream part of the Red River delta has been set-up to study tides with the 2D spectral TUGO model and basin-scale circulation with the 3D Symphonie model. The PhD thesis of C. Tran (ongoing) at LOG (Wimereux, France) aimed at better understanding the circulation in the Gulf of Tonkin from tidal to annual scale using the High Frequency (HF) radars installed near the coast.

This section will start with the definition of what a river plume is and the scientific and practical interests of their study. It continues with a description of the regional context of the Gulf of Tonkin and Red River delta, as well as an account of previous studies in the area. Finally, we announce what will be our objectives and approaches in this work, and the outline of this manuscript.

I.2 River plume

I.2.1 What is a river plume?

A river plume can be defined in a general way as the region of the coastal ocean where its properties and dynamics are affected by the river runoff (Horner-Devine et al. 2015). Though the river runoff is small compared to the whole ocean water volume, it can impact both the physics and biogeochemistry of the coastal ocean depending on the discharge, the properties of the ocean area (bathymetry, bottom roughness) and external forcing (air-sea fluxes, open ocean influence). (2015).

The dynamics of a prototypical plume has been studied by several authors [START_REF] Fong | The alongshore transport of freshwater in a surface-trapped river plume*[END_REF][START_REF] Hetland | Relating river plume structure to vertical mixing[END_REF][START_REF] Horner-Devine | Laboratory experiments simulating a coastal river inflow[END_REF]. In a simplified coastal numerical simulation, with relatively high discharge and narrow river mouth, [START_REF] Hetland | Relating river plume structure to vertical mixing[END_REF][START_REF] Hetland | The effects of mixing and spreading on density in near-field river plumes[END_REF] The source water is the area before the liftoff (the point where the buoyant outflow detaches from the bottom). The near field starts from the liftoff and forms the buoyant layer. In this area, the water flows like a jet. Therefore, the limit of this area is where the supercritical flow ends. Then, in the mid-field, the plume is affected by both the momentum of the river flow and the effect of Earth rotation. Under low wind conditions, it can form a recirculating bulge region. Finally, when the plume loses the momentum of the initial runoff, and is driven mostly by buoyancy, earth's rotation, wind and ambient currents, but is still distinguished with the ambient water, it is the far field plume.

Depending on the runoff, the characteristics of the river mouths, combined with the atmospheric and oceanic conditions, the structures of the river plume are different. Not every single river plume has all the features of the prototypical plume as described above. Garvine (1995) suggests a classification of the river plume, based on Kelvin number, equal to the ratio of the width of the river mouth to the Rossby radius of deformation. Plumes with high Kelvin number (>>1) are considered large scale plumes, while small scale plumes are characterized by small (<<1) Kelvin numbers. [START_REF] Huq | The role of kelvin number on bulge formation from estuarine buoyant outflows[END_REF] shows that the bulge is formed if the estuary is narrow compared to the Rossby radius (i.e Kelvin number < 1).

Horner-Devine et al. ( 2015) divide the plume types into 6 general categories. The prototypical plume (A), as described above, is characterized by a large runoff, narrow river mouth, and its dynamics is affected sharply by Earth's rotation. The non rotational plume (B), in contrast, is less impacted by Earth's rotation due to small runoff or the river mouth is located near the Equator. Wide estuary plume (C), as its name implies, has a wide mouth and small runoff. Inside the river, the runoff flows in one side of the river due to Earth's rotation. The fourth type, angled inflow plume (D), has a river with a small angle to the coast. In delta plume (E), the runoff flows through several river mouths and has a joint plume from different mouths. Finally, region of freshwater influence (ROFI), appears in shallow marginal seas ; the plume has a strong interaction with the bottom.

I.2.2 Why do we need to study the river plume?

Rivers are vectors of exchange between the continent and the ocean. They bring freshwater that impacts the buoyancy of the coastal waters and therefore the coastal dynamics; the rivers impact is sometimes very local, sometimes at the basin scale (as for example in the case of the Amazon). Rivers carry matter and sediments that influence the coastal morphology (bathymetry, coastline). They also carry chemical species either in dissolved or particulate forms (for instance in sediments); some of them are natural while others have an anthropogenic origin from agriculture and industrial activities. Both influence the water quality and the ecosystems; for instance loadings of nutrients (nitrogen, phosphorus …) mainly stemming from watersheds with intensive farming can lead to algal bloom. At last, river flows bring along other kind of contaminants such as plastic debris.

Therefore, there is a need to better understand the fate of the river water from the estuaries to the ocean, especially in the mid-field and far field. Using a fuzzy c-mean clustering, [START_REF] Delaval | A fuzzy classification of the hydrodynamic forcings of the Rhone River plume: An application in case of accidental release of radionuclides[END_REF] classify wind and Rhone River discharge into 6 scenarios. From that, simulations have been done in each scenario to identify the plume pattern in different forcing conditions to predict the main Rhone River plume patterns and consequences in case of radionuclides accidental releases. Using numerical simulation, Menesguen et al. (2018) identify the source of eutrophication in the Bay of Biscay and the English Channel and design a strategy of nutrient load reduction within European institutions directives framework.

As shown above, understanding the variability of the river plume in different scales is the first step toward the study of the dispersion of the possible contaminations and toward the design of strategies for monitoring and managing the water quality and the health of ecosystems. These are particularly crucial issues in densely populated areas, such as many deltaic regions of Southeast Asia, including the Gulf of Tonkin.

I.3 The Gulf of Tonkin and the Red River

The

I.3.1 Atmospheric conditions

The wind conditions in GOT depend on the monsoon. In the winter (October to February), the wind direction is from the northeast, with the velocity from 4 -8m/s (Figure I.3). The wind is strongest in December, with the maximum value of 8m/s in the center of the gulf. The summer monsoon is shorter, from June to August, with less intense wind (2 -6m/s). Between the two monsoons is the transition period with easterly wind. The wind is weakest in September, with velocity less than 2 m/s over the whole gulf. Due to the monsoon, the river plume undergoes a seasonal variation which will be presented in Chapter 4. The monthly mean precipitation in the GOT is shown in Figure I.4. From January to May, it is less than 5mm/day, with higher values located near the Vietnamese coast. Then, from June to August, the precipitation pattern changes with the highest values located in the north of the gulf. Over most of the area, the precipitation ranges from 5 -10mm/day, with a peak of around 15mm/day in July. In September, the pattern changes again with the maximum value located near the Vietnamese coast. The precipitation is mostly higher than 10mm/day, with the local maximum of more than 15 mm/day from 18 to 20N. The high precipitation from July to September is expected to significantly impact the sea surface salinity, with a possible ambiguity when trying to identify the low salinity water of the river plume as discussed in Chapter 4. From October, the precipitation decreases again. In December, it is less than 5 mm/day with a local maximum of 8 mm/day near the coast between 17 -19N.

I.3.2 Circulation

The ocean circulation of the Gulf of Tonkin (hereafter GOT) has been examined by several authors (e.g. Ding et al., 2013;Gao et al., 2014;Rogowski et al., 2019; for a recent review see Piton et al., 2021). Using model results and observational data, they all agree that the winter circulation at the gulf scale is cyclonic and driven by wind. However, in summer, the circulation is still a hot debate. Based on model results and 10-m current measurements, Ding et al. (2013) claimed that the flow in the northern Gulf of Tonkin features a cyclonic circulation. By conducting different numerical experiments, they also showed that this cyclonic circulation is induced by tide and can be accelerated by the wind. Gao et al. (2013) showed that the current in the southern gulf is anti-cyclonic and driven by the general South China Sea circulation.

Using a 3D numerical model, HF radars and in situ data, Rogowski et al. (2019) showed that the current pattern is complex in the north of the gulf. Meanwhile, there is a strong surface water flow from the gulf to the South of Hainan island in June and strong Southward Vietnamese coastal flow in September in the south. From high-resolution model and pluriannual simulations, Piton et al. (2021) found a basin-scale anticyclonic circulation in summer. They also suggest that the surface circulation is mainly ageostrophic, as a consequence of the monsoon wind forcing, except along the Vietnamese coast where the southward coastal current has a dominant density-driven component.

The flow through Hainan Strait during the summer monsoon is also discussed in several papers but is still an open question. Shi et al. (2002) and [START_REF] Chen | Impact of the water input from the eastern Qiongzhou Strait to the Beibu Gulf on Guangxi coastal circulation[END_REF] share the same opinion that the transport through Hainan Strait is inflow (or westward or flow into the Gulf) while Gao et al. (2013) showed that the mean current flows out of the gulf in the late spring and summer and into the gulf in other seasons. They also noted that the daily variations of the current are very strong so that the daily flow direction could vary independently of the season. [START_REF] Zhu | Tidal and residual currents in the Qiongzhou Strait estimated from shipboard ADCP data using a modified tidal harmonic analysis method[END_REF] showed that the residual current is inflow in spring 2013 and controlled by tidal rectification.

The link between the flow direction in Hainan Strait and the circulation in GOT is studied by Wu et al. (2008). They showed that in the summer, if the flow through Hainan Strait is inflow, the circulation is mainly cyclonic and vice versa.

I.3.3 Tides

Tides in the South China Sea in general and GOT in particular have been studied for many years. [START_REF] Wyrtki | Scientific Results of marine Investigations of the South China Sea and the Gulf of Thailand 1959-1961[END_REF] shows that 4 main constituents (K1, O1, M2, S2) are enough to explain the tidal signal in the South China Sea. Until recently, there were only a few studies about the tidal characteristics in the South China Sea in general [START_REF] Cai | Tide model evaluation under different conditions[END_REF][START_REF] Fang | Numerical simulation of principal tidal constituents in the south china sea, gulf of tonkin and gulf of thailand[END_REF]) and in the Hainan Strait [START_REF] Chen | Numerical study of the tides and residual currents in the Qiongzhou Strait[END_REF]. More recently, using a high resolution model with a horizontal resolution of 4.5 km and satellite and in situ observations, the tidal dynamics in the GOT has been examined by Nguyen et al. (2014). The result shows that the diurnal constituents (K1 and O1) are predominant in the GOT. The largest amplitudes are located in the north of the gulf, with amplitudes of more than 90 and 100 cm, respectively (Figure I.5). South of 18N, the amplitude decreases to under 30 cm for both diurnal constituents. M2 and S2 have similar patterns, with higher valuse in the west, along the coast of Vietnam, and the northeast coast of China. While the maximum value of M2 is 70 cm, the maximum value of S2 is 20 cm. The Red River (RR) is the second-largest river system in Vietnam after the Mekong, and is the main river that feeds the gulf (Figure I.6). It originates in the Yunnan Chinese province and owes its name to the color of its waters which carry iron-rich sediments [START_REF] Vinh | Suspended sediment dynamics in Red River distributaries and along the Red River delta: focus on estuarine processes and recent balances[END_REF]. The RR system is formed by 3 tributaries: Da, Thao and Lo Rivers which join upstream the city of Son Tay (about 150 km from the sea). A few kilometers upstream of Ha Noi, they split into the Red River and the Thai Binh river, which in turn split into a network of connected distributaries: the Red River delta. A complex system of channels has also been built in the delta for irrigation and drainage purposes.Daily RR runoff, obtained from the National Hydro-Meteorological Service (NHMS) of Vietnam at Son Tay hydrological station, shows that the mean runoff from 2008 -2018 is approximately 3000 m3/s, accounting for more than 60% of the total discharge to the gulf. The runoff undergoes a strong seasonal cycle with high discharge during the rainy season (from June to September) and low discharges in the dry season (October to May). Interannual variability affects the hydrological regime as well, but regulation intervenes to redistribute water in the channels and tributaries (see for instance [START_REF] Luu | Hydrological regime and water budget of the Red River Delta (Northern Vietnam)[END_REF].

The sediment transport from the Red River varies greatly at seasonal and interannual time scales. According to [START_REF] Lu | Sediment budget as affected by construction of a sequence of dams in the lower Red River, Viet Nam[END_REF], the sediment load at Son Tay in [START_REF] Pierson | Relationship between the attenuation of downwelling irradiance at 490 nm with the attenuation of PAR (400 nm-700 nm) in the Baltic Sea[END_REF], 2009, and 2010are 30, 12.8, and 10 Mt/yr, respectively. Wei et al. (2021) show that the mean sediment load from 2008 -2013 is 11.6 Mt/yr.

The RR delta (RRD) is formed from 11 provinces (Prime Minister, 2013), which corresponds to a total population of 22.5 millions (Vietnam General Statistic Office. 2020) and an area of 21 200 km2. The average population density of the RRD is around 1060 inhabitants/km2, which is nearly four times higher than the national average for Vietnam. The very high population density exerts a strong pressure on the rivers and their environment. For instance, large contaminant spills are expected, due to farming and industrial activities.

One of the pollution problems is the emergence of microplastics in the ocean. Microplastics are plastic particles ranging from 1 to 5000 μm long [START_REF] Frias | Microplastics: Finding a consensus on the definition[END_REF]. The measurement made by [START_REF] Strady | Baseline assessment of microplastic concentrations in marine and freshwater environments of a developing Southeast Asian country, Viet Nam[END_REF] showed that the microplastics concentration in Red River is around 2.3 items/m3, while in the ocean, it is 0.4, 0.7 and 0.8 items/m3 in Cua Luc, Do Son and Ha Long Bays, respectively. Better understanding the variability of the river plume can help us to identify the pathway of the microplastics from river to ocean.

South of the RR, there are 3 more rivers: Ma, Yen, Lam. The summed discharges of these rivers and RR can account for around 90% of the total runoff in the GOT.

I.3.5 Red River Plume

In spite of its importance, to date, studies focusing on the RR plume are still scarce.

van [START_REF] Van Maren | Water and sediment dynamics in the Red River mouth and adjacent coastal zone[END_REF] studies the water and sediment dynamics inside the Ba Lat mouth and its near field plume. The measurements at 5km inside the river mouth and at 2 km seaward, performed in the low and high discharge period in 2000 show that the temperature of the fresh river plume was usually 1-3° lower than ambient seawater. The freshwater plume (with salinity less than 6) is separated from the ambient water at 3km from the river mouth. He also suggests that in the dry season, the plume is subcritical, characterised by strong coastal current, while in the wet season, the plume is supercritical with large offshore spreading.

Using 2-years numerical simulation, Gao et al. (2013) showed that the coastal plume was found near the northern and western coasts of the gulf in winter while spreading eastward and offshore in summer. They suggest that the difference between the 2 seasons due to the monsoon wind conditions, which is downwelling favorable in the winter and upwelling favorable in the summer. The 10 -years simulation of Rogowski et al. (2019) show similar results : the southwesterly wind direction which is prominent during the summer monsoon is the main mechanism preventing the downcoast advection of the RR plume.

Using the CTD data located from 2 -25 km from the coast in 2 years 2015-2016, Rogowski et al. (2019) and Tran et al. (2021) tried to examine the vertical structure of the RR plume from the stratification of salinity in different time of the year. In September, the plume thickness is 5m and extend beyond the survey line (25km from the coast). In the two following months, the plume is deeper and detaches from the coast. In January, the plume is well mixed, due to the cooling of the coastal water. It is also the time when the runoff is lowest. The low stratification appears until March. In April, the thin plume is observed on the surface. There is no plume signature in June and July, suggesting that the river plume is advected offshore under southwesterly wind conditions. To my knowledge, the existing studies still have some limitations to describe and explain the overall picture of the RR plume variability. In the study of van Maren, D. S. ( 2007), the role of wind is not examined because the wind velocity is very low during the measurement. The simulation of Gao et al. (2013) is run with climatological mean runoff, therefore cannot show the interannual variability of the plume. The simulation of Rogowski et al. (2019) has aspatial resolution of 6.5km, which is likely insufficient to resolve the dynamics associated with the RR plume.

I.4 Objectives and approaches

In this study, the main objective is to provide a better understanding of the variability of the RR plume in the GOT in different scales. The main scientific questions which we wish to address are the following:

1. How can high-resolution numerical simulations reproduce the Gulf of Tonkin dynamics and the Red River plume variability at different scales in the study years? 2. Which criteria can be used to identify the RR plume in the GOT from the model outputs? 3. How can we formally characterize the sequence of regimes which the RR plume goes through in the study years? Can we link those regimes to specific processes and forcing functions? 4. How can we assess the sensitivity of model variables and in particular on the variability of the river plume, to the main sources of uncertainty, in particular the atmospheric forcing uncertainties?

To address those questions, we will follow an approach involving several tools:

1. We will set up and run an upgraded version of the model configuration created as part of the thesis of V. Piton. 2. The river plume will be identified using different methods, allowing us to select the most appropriate one given our purposes. 3. The RR plume regimes will be identified and linked to physical processes with the help of K-means clustering analysis. 4. We will run an ensemble of simulations assessing the model response due to perturbations of the wind forcing. We will analyze the variability of the resulting uncertainties over space and time of the ocean variables. The K-means clustering approach will be used on ensemble members, giving access to an assessment of the robustness of the results of 3. above.

I.5 Outline of the thesis

The structure of the thesis is as follows. Chapter 1 presents an overview of the river plume and some past studies of the river plume in the different regions. Then, the physical characteristics and an overview of the past studies in GOT is presented. In Chapter 2, the SYMPHONIE model, model configurations, data sources and the unsupervised learning algorithm used to classify the main pattern of the plume (K-means) are described. The model is then evaluated against several observational data sets, including HF radar current data from the results of C. Tran thesis, in Chapter 3. The impact on the dynamics of a new parameterization of the light attenuation in the water column is also presented.

Chapter 4 contains an article published as a part of this thesis. To identify the plume in the GOT, several methods are compared which allows us to select the most appropriate one given our purposes. On that basis, the river plume variability at different time scales is examined with the help of K-means clustering analysis, illustrating the effect of key physical processes. Then, a discussion on the plume classification and a description of the receiving basins from the different river systems in the GOT are presented. The end of this chapter provides some further discussions on the model boundary conditions and on the sensitivity of the clustering algorithm.

Chapter 5 aims at assessing the robustness of the model in general, and of the river plume estimate of Chapter 4 in particular, to the uncertainty of the wind forcing. To do that, an ensemble simulation approach is used. Firstly, the uncertainty of the wind forcing is assessed by comparing it to satellite data. A pseudorandom combination of bivariate wind EOFs is then used to perturb the wind field within that uncertainty. The sensitivity of the main model variables to the perturbations is then assessed by calculating the ensemble spread and the distributions of the variables of interest. Then, the sensitivity of the river plume to the stochastic forcing is presented. To assess the effectiveness of the ensemble configuration at reproducing the error, the ensemble is verified using the available observational datasets.

Finally, Chapter 6 summarizes and gives overall conclusions and future perspectives.

CHAPTER II: METHODOLOGIES AND DATA

*** Most of the investigations, tools and approaches in this Ph.D. thesis use numerical modeling. In this chapter, I will first describe the principle of the SYMPHONIE model. Then, the configuration of the study area is presented. Finally, a list of the available observational data which are used for model validation is provided.

II.1 SYMPHONIE model

SYMPHONIE is a numerical model that solves the primitive, Boussinesq, hydrostatic equations of the ocean circulation (Marsaleix et al. 2006(Marsaleix et al. , 2008) ) on a curvilinear bipolar (Bentsen et al. 1999) Arakawa C-grid with sigma or hybrid vertical levels. SYMPHONIE has been developed to model the coastal circulation from regional scales to shelf seas and nearshore areas. In particular, SYMPHONIE allows to accurately represent tides, intertidal areas, and the oceanic response to high-frequency atmospheric forcing. It can run on highly variable polar or bipolar grids in order to refine the resolution locally, for instance at river mouths or inside estuaries. It has been extensively used in the Mediterranean Sea to study a wide range of physical processes, from the large scale circulation, to the deep convection in the western basin, dense water cascading, river plume, impact of storms, wave-induced currents, sediment transport and impact of the circulation on the biogeochemistry (see https://sirocco.obs-mip.fr). It has also been used in the Bay of Biscay and Gironde estuary, Gulf of Guinea, northern Pacific ocean and southeast asian seas. Since the PhD of V. Piton, it has been implemented in the Gulf of Tonkin. This section is devoted to the description of the model equations and numerical implementation. where 𝐼 $ is the part of the solar radiative forcing that penetrates to depth z, 𝐶 % is the water heat capacity. The horizontal diffusion fluxes for temperature and salinity are not explicitly simulated; instead, the diffusivity of the numerical scheme for the advection is used (QUICKEST scheme). The eddy diffusivity coefficient is determined by the turbulence closure scheme.

II.1.1 Governing equations

II.1.2 Turbulence closure scheme

Coefficients of eddy diffusivity (𝐾 & ) and viscosity (𝐾 # ) are computed using turbulence closure schemes. In this study, the K-𝜀 turbulence closure scheme is used, with the implementation described in Michaud et al. (2012).

The vertical eddy diffusivity 𝐾 ! and viscosity 𝐾 " are defined as: .% ) can be expressed, following [START_REF] Craig | Modeling wave-enhanced turbulence in the ocean surface layer[END_REF], as:

𝐾 ! = #2𝐸 # 𝑙 # 𝑆 ! Eq.
𝐹 % = 100 K 𝜏 ) 𝜌 L , "
Eq. 2.15

where 𝜏 $ is the wind stress computed by the bulk formula using the ECMWF variables (wind velocity, temperature and relative humidity).

The surface current condition affected by the wind stress τs is determined as follows:

𝜌 ! 𝐾 # K 𝜕𝑢 𝜕𝑧 , 𝜕𝑣 𝜕𝑧 L = Q𝜏 )0 , 𝜏 )1 R Eq. 2.16
The pressure is forced by the atmospheric pressure at the surface provided by ECMWF.

The heat flux at the ocean-atmosphere boundary is computed from 4 components.

𝐾 # 𝜕𝑇 𝜕𝑧 = 𝑄 2 + 𝑄 3 + 𝑄 / + 𝑄 4
Eq. 2.17

With:

• 𝑄 ( : solar radiation

𝐾 # 𝜕𝑆 𝜕𝑧 = (𝐸 2 + 𝑃𝑟 2 )𝑆 Eq. 2.18
with 𝐸 ( is the evaporation rate and 𝑃𝑟 ( is the precipitation rate.

The solar radiation, the downwelling longwave radiation and the precipitation are provided by ECMWF analyses. The longwave radiation from the surface, turbulent heat fluxes and evaporation are computed using atmospheric variables from ECMWF analyses, the sea surface temperature simulated by SYMPHONIE and the bulk formulae of Large and Yeager (2004).

II.1.3.2 Bottom boundary

At the bottom, the boundary condition approximation for 𝐸 ' corresponds to the equilibrium between energy production and dissipation (Michaud et al., 2012) The open-boundary conditions are described in Marsaleix et al. (2006) and Toublanc et al. (2018, their appendix). At an open ocean boundary, the equations write differently for the barotropic and baroclinic velocities. Here, the barotropic velocities are defined as the depthaveraged velocities, while the baroclinic ones are the anomalies with respect to the vertical mean. For simplicity, we suppose that the domain of integration is 1D (along x) with open boundaries at x=0 and x=L:

For the barotropic current (𝑢 :, 𝑣̅ ):, the Flather type condition is applied to the surface elevation η:

𝜂 -𝜂 : = -] 𝐻 𝑔 Q𝑢 : -𝑢 : :::R

at x = 0 Eq. 2.23a
𝜂 -𝜂 : = ] 𝐻 𝑔 Q𝑢 : -𝑢 : :::R

at x = L Eq. 2.23b
where the f index indicates the forcing (large scale circulation and tidal circulation).

The tangential transport at the boundaries is computed using the Neumann condition:

𝜕𝐻Q𝑣̅ -𝑣 : :::R 𝜕𝑥 = 0 at x = 0 and x = L Eq. 2.24

The transport component normal to the boundary is deduced from the continuity equation, from η (Eq. 2.23) and from the tangential transport (Eq. 2.24).

For the baroclinic currents (u, v), a radiative condition is applied on the difference between the SYMPHONIE variables (Ф) and the forcing variables (Фf): Both the barotropic and baroclinic velocities are relaxed towards the external forcing within a sponge-layer in order to reduce possible reflection of the outgoing flow. In the sponge layer, the relaxation time scale decreases from a fixed value at the boundary to zero at the most internal point of the sponge-layer inside the integration domain.

For the scalar variables (T, S), the open boundary conditions use a hybrid scheme depending on the flow direction. If the flow is incoming, the temperature and salinity (T, S) are deduced from the forcing (Tf, Sf) of a larger scale model; a relaxation to (Tf, Sf) is also applied in the sponge layer with a time scale depending on the forcing current amplitude. In the case of outflow, these variables are issued from the simulated domain.

II.1.3.4 Boundary condition at the river input point

At the river input point, the river boundary current is consistent with the runoff of the river as in Reffray et al (2004), while the vertical shape of the current is given by a gradient condition allowing the current to adjust by itself to the dynamics of the river downstream (rather than arbitrarily prescribing a possibly inappropriate vertical profile). The boundary condition is applied in two steps.

Step 1 is given by the above-mentioned gradient condition, leading to a provisional solution of the current:

.= * .0 = 0 if current in downstream direction of the river u* = 0 otherwise Eq. 2.26
where x is the direction of the river flow and u* a provisional solution of the current.

Step 2 ensures that the total flow is consistent with the expected runoff from the river. In practice, the river boundary current is proportional to u*, i.e. 𝑢 = 𝛼𝑢 * where 𝛼 is a constant used to satisfy the constraint on the prescribed river runoff:

𝛼𝛥𝑦 [ 𝑢(𝑧) * 𝜂 -ℎ 𝑑𝑧 = 𝐹 Eq. 2.27
where 𝛥𝑦 is the cellbox width at the river input point and F the river runoff.

II.1.4 Tides

The tidal signal is made of different constituents corresponding to different frequencies. The tidal surface elevation and barotropic currents can be written as follows (see for instance Pairaud et al., 2008): where η0, u0, v0 and Gη, Gu, Gv are the amplitudes and phase lags for sea surface elevations and currents of the tidal constituent k at the longitude λ and the latitude φ. The nodal factors fk and uk are considered as constants [START_REF] Doodson | The analysis of tidal observations[END_REF], ω is the tidal frequency and V0 is a constant related to the reference time t0.

𝜂(𝜆, 𝜑, 𝑡) = ] 𝑘 𝑓 𝑘 𝜂 0,𝑘
In coastal models, tides are prescribed at the open boundaries, usually as harmonic constituents (amplitude and phase) for the main frequencies in the areas of study (the tidal signal is then reconstructed at each time step using Eq. 25). The prescribed fields come from global atlases (such as FES2014) or dedicated simulations (as discussed for instance in Toublanc et al., 2018).

In SYMPHONIE, the tidal forcing also consists of astronomical and loading potentials (respectively ΠA and ΠL) (see Pairaud et al., 2008, their appendix). The tidal potential due to astronomical effects ΠA [START_REF] Hendershott | The Effects of Solid Earth Deformation on Global Ocean Tides[END_REF] is taken into account as follows:

𝛱 > = (1 + 𝑘 " + ℎ " )𝑓𝑎 × f𝑣 ! 1 -3𝑠𝑖𝑛 " (𝜑) 2 + 𝑣 + 𝑠𝑖𝑛(2𝜑) + 𝑣 " 𝑐𝑜𝑠 " (𝜑)l × 𝑐𝑜𝑠Q𝜔(𝑡 -𝑡 ! )R + 𝜈𝜆 + 𝑉 ! + 𝑢
Eq. 2.29 with:

𝑘 ! , ℎ ! : Love numbers 𝑎: the equilibrium amplitude of the tidal constituents considered, v0, v1, v2, and ν: parameters depending on the type of tide. In case of long period tide v0 = 1 and ν = 0, of diurnal tide v1 = 1 and ν = 1, of semi diurnal tide v2 = 1 and ν = 2.

The loading potential ΠL (LSA potential) is expressed as:

𝛱 3 = 𝑓 𝛱 ! (𝜆, 𝜑)𝑐𝑜𝑠(𝜔(𝑡 -𝑡 ! ) + 𝐺 ? (𝜆, 𝜑) + 𝑉 ! + 𝑢) Eq. 2.30
The hydrostatic pressure force in the horizontal momentum equation is determined by:

- 𝛻𝑝 𝜌 ! = - 𝑔 𝜌 ! 𝛻 7 (𝜌 -𝜌 ! )𝑑𝑧′ $ % -𝑔𝛻(𝜂 -𝛱 > -𝛱 3 )
Eq. 2.31 A wetting and drying scheme is used to represent intertidal areas.

The model computes an online harmonic analysis on the sea surface elevation and depthaveraged currents. A first run is performed to compute the tides. Then the result of the harmonic analysis from this first run is used to detide the daily averages of the sea surface elevation and horizontal current that we use in our studies.

II.1.5 Discretization

II.1.5.1. Spatial discretization
The SYMPHONIE model equations are solved using the finite difference method (Marsaleix et al., 2008) on the Arakawa staggered C-grid [START_REF] Arakawa | Vertical differencing of the primitive equations in sigma coordinates[END_REF]) (Figure II.1).

Horizontal velocity components (u, v) are defined at the center of the edges of the grid, whereas the scalar components (temperature, salinity, surface elevation) are defined at the center of the grid cell. The vertical velocity and kinetic energy are calculated at the center of the grid edges at each vertical level. [START_REF] Piton | Du Fleuve Rouge au Golfe du Tonkin: dynamique et transport sédimentaire le long du continuum estuaire-zone côtière[END_REF]. The i, j, k indices correspond to x, y, z direction, respectively. φ represents the scalar fields. Turbulent fields are defined at the same positions as the vertical component of velocity.

For the vertical grid, SYMPHONIE provides several options. In this study, due to the shallow bathymetry in the GOT, the terrain-following 'sigma' coordinates system is used. It allows us to model with high vertical resolution the coastal shallow area where the plume spreads while saving the computing resources in the deep area.

The vertical advection is computed using a centered scheme. Horizontal advection and diffusion of temperature and salinity are computed by the QUICKEST scheme [START_REF] Leonard | A stable and accurate convection modelling procedure based on quadratic upstream interpolation[END_REF]. Horizontal advection and diffusion of momentum are each computed with a fourth order centered biharmonic scheme as in Damien et al (2017). The biharmonic viscosity of momentum is calculated according to a Smagorinsky-like formulation derived from Griffies and Hallberg (2000).

II.1.5.2. Temporal discretization

The forward-backward (F-B) time-stepping scheme is used for the temporal discretization of the equations. In this explicit time-stepping methods, the variable value at time t+1 is deduced by the variable value at time t:

𝐹′ = 𝐹 @ + 𝛥𝑡𝑓(𝐹 @ )
Eq. 2.32a

𝐹 @A+ = 𝐹 @ + 𝛥𝑡𝑓(𝐹′)

Eq. 2.32b

The F-B time-stepping application and chronology of variables in SYMPHONIE are described in more detail in [START_REF] Marsaleix | 3D phase-resolved wave modelling with a nonhydrostatic ocean circulation model[END_REF].

The model uses a time-splitting method with a double loop on time: the model solves first the external mode (SSH and barotropic horizontal velocities) then solves the internal mode (baroclinic velocities, tracers and SSH). The external mode is solved with a very small time step consistent with fast gravity waves (such as tides).

II.2 Model configuration

In this study, the general configuration has been initially set-up for the PhD of V.

II.2.1 General configuration

The model domain covers the GOT area, with 429 * 600 grid points. Due to the variable horizontal grid, the coastal area near the Red River (RR) mouths has a fine horizontal resolution of 300m, while near the open boundary, the grid size can increase to 4500m (Figure II.2). The vertical discretization consists of 20 sigma levels. As the coastal area adjacent to the RR delta is characterized by a complex topography, with many islands and islets, a considerable effort has been devoted to the construction of the model bathymetry and shorelines and is described in Piton et al. (2020). In particular, the bathymetry is reconstructed from GEBCO 2014 combined with other sources and field surveys. The bottom drag coefficient follows a logarithmic law (Blumberg and Mellor, 1987), depending upon the bottom roughness length (z0). Using the ROMS model, Minh et al. (2014) showed that the tidal error for O1 and K1 is lowest with z0 around 1 mm. Piton et al. (2020), in their sensitivity study using a 2D configuration of the T-UGOm model, shows that the sum of errors for the 4 main tidal constituents (O1, K1, M2, S2) is lowest with constant z0 equal to 0.9 mm. Furthermore, the error can be even more reduced with spatially varying z0. However, these results do not hold in the 3D SYMPHONIE configuration: in Piton (2019), she confirmed that using z0 = 1 mm gives a better solution in simulating the tide in. Therefore, in this study, z0 is set up equal to 1 mm. At the surface, boundary conditions are provided by an atmospheric model and fluxes of momentum, heat and freshwater are computed using the bulk formulae of Large and Yeager (2004). Operational ECMWF analyses (with a spatial resolution of 1/8°) are used to provide 3hourly wind, precipitation, solar energy, atmospheric temperature, dew-point temperature, surface pressure. The land points of the atmospheric variables are not masked during the interpolation processes.

As the present study focuses on the fate of continental water into the coastal ocean, a specific effort has been deployed on the river runoffs implementation; section II. The daily RR runoff is obtained from the National Hydro-Meteorological Service (NHMS) of Vietnam at Son Tay hydrological station which is located at the apex of the RR delta. For the southern rivers, since the daily runoffs were not available, the monthly climatological runoffs from NHMS is prescribed. The monthly runoffs of the northern rivers are obtained from the data given by Gao et al. (2013).

Figure II.4 shows the discharge from 3 river systems. In general, the discharge from the RR is the largest. On average, the runoff of the Red River system alone accounts for 60% of the total runoff to the gulf. In detail, the average discharge rate for the Red River system in low (December, January, February) and high (July, August, September) discharge periods equal 1632 m 3 /s and 4959 m 3 /s, respectively. It also shows interannual variations. While the peak in 2012 is approximately 15 000 m 3 /s, the peak in 2011 is just over 6 000 m 3 /s.

For the SR, runoff in high and low discharge periods is 365 and 2043 m 3 /s, respectively. For the NR, they are 164 m 3 /s and 1103 m 3 /s. Since these values are monthly averaged climatological data, it does not have interannual variability and can be a source of uncertainty when assessing the plume area.

The ratio between high and low discharge seasons for RR is 3, while it is 5.6 for the SR and 6.7 for the NR. The lower ratio for the Red River may be explained by water flow regulation by several hydrological dams upstream. On average, the total runoff from RR and SR accounts for approximately 90% of the total runoff in the GOT.

II.2.2.2 River configuration

In most modeling studies to date, the RR was set up with only one mouth (or input grid point) using monthly or annual climatological runoff because more realistic runoff at the hydrological stations was unavailable (Ding et al. 2013, Gao et al. 2013). In this study and in Piton et al. (2021), the river condition is configured as realistically as possible.

Using the river discharge at the upstream and the water level at the river mouths, Vinh et al. (2014) attempted to model the Red River network and calculated the distribution of water discharge from Son Tay hydrological station into different river mouths. The results showed that the runoff is not equi-distributed. In this model configuration, the discharge is distributed across the seven distributaries based on their results: Bach Dang (7%), Cam (13.2%), Van Uc (14.5%), Thai Binh (6.4%), Ba Lat (30.3%), Ninh Co (5.6%), Day (23.0%). These 7 mouths lie along 150km of the coastline, and are expected to represent a better distribution of the plume than a single input point.

Furthermore, each river mouth is connected to a channel in the model to allow the water to enter the coastal ocean with realistic salt and temperature properties and realistic stratification. The length of the channels is chosen to exceed the saltwater intrusion, which is approximately 30 km from the mouth [START_REF] Nguyen Thi Hien | Upstream effects on salinity dynamics in the Red River Delta[END_REF]. In detail, the length for each channel is: Bach Dang (57 km), Cam (78 km), Van Uc (221km), Thai Binh (115km), Ba Lat (120km), Ninh Co (72km), Day (72km). At the most upstream point of the channel (input point), the salinity of the river flow is set to 0. The temperature varies seasonally from 17°C (in February) to 29°C (in August) to be consistent with observational data from a river in the north of Vietnam.

In the model, the channel can be configured automatically using some parameters (basic mode) or can be defined in each grid point (advanced mode). In the basic mode, the channel is configurable by the grid length (L), grid width (W) and the height (H) at the most upstream point (Figure II.5).

The entrance of the channel (downstream) is connected with a grid point on the ocean. The elevation of the river bed with respect to the reference level (z=0) at each grid point along the channel is linearly interpolated from the depth of the downstream point and H. We do not define the depth of the channel, but we assume that there is no overflow. It means that the water height only depends on the discharge and the width of the channel. Similarly, the width of each grid point along the channel is interpolated from W and the width of the ocean point. The length is configured similarly. The channel is straight, and it does not represent the shape of the distributary. In the simulation without tide, the salinity intrusion seems more stable within 1 day. The salinity in both time frames shows strong stratification at around 25km from the river mouth. It also shows the changes of SSH over the day in GOT_REF, while in GOT_NOTIDE the change in SSH is much less significant. The important thing is that in both cases, with and without tide, the salinity at the river mouth is not zero and the implementation of the channel can help create better salinity mixing at the river mouth.

Passive tracer configuration

There are several ways to identify the plume. It can be detected using the difference between the salinity of the riverine water and the ambient water (e.g. Burla et al. (2010), [START_REF] Nguyen-Duy | A model study of tide-and wind-induced mixing in the Columbia river estuary and plume[END_REF] for Columbia river, Falcieri et al. (2014) for Po river ), or using the satellite data thanks to the difference in temperature, turbidity, salinity, of the riverine water (e.g. [START_REF] Brando | High-resolution satellite turbidity and sea surface temperature observations of river plume interactions during a significant flood event[END_REF][START_REF] Hopkins | Detection and variability of the Congo River plume from satellite derived sea surface temperature, salinity, ocean colour and sea level[END_REF].

In GOT, the precipitation has seasonal variations, as shown in Figure I.4, with large precipitation along the coast during the high discharge period. Therefore, using surface salinity to identify the river plume could be biased by the freshening due to precipitation. In order to simulate the pathways of the river water into the GOT, in this study, passive tracers are used. Tracers act as dyes, i.e., they do not affect the dynamics. This approach has been used in several studies. [START_REF] Jia | Summertime Connecticut river water pathways and wind impacts[END_REF] use a tracer to simulate the pathway of the Connecticut River in the summer under different atmospheric conditions. The distribution of tracer in different areas allows to identify the progression of the river water at the tidal time scale. [START_REF] Gong | Plume-to-plume interactions in the Pearl River Delta in winter[END_REF] use several tracers to examine the plume-to-plume interaction in the Pearl River plume in the winter, when the runoff is highest. The results show that different plumes can affect each other and the existence of a plume can limit the transport of others. Wang et al. (2014) use a passive tracer to study the variability of the Copper River plume by calculating the evolution of the plume volume over time and plume occurrence frequency in different periods. This is the first time that tracers are used for a pure dynamical study with SYMPHONIE. In total, there are 3 tracers (or three colors) meant to distinguish the inputs from the different river systems. The first tracer is added to the runoff of all the Red River distributaries at constant concentration (100 arbitrary unit/m 3 ). The second one is added at the mouths of the rivers in the south and the third one at the mouths of the rivers in the north of RR, with the same concentration (100 arbitrary unit/m 3 ). Since tracers are injected with the runoff, they are also present in the 1 year of spin-up. The boundary conditions for the passive tracer (which is put directly into the river water) are applied in the same way as the salinity, except that the incoming concentration is always zero at the lateral boundary. Therefore, an appropriate heat flux can improve the model result.

II.2.3 Light attenuation parameterization

One of the parameters controlling heat fluxes is the shortwave attenuation depth (hereafter SWAD or L2). It characterizes the ability of the shortwave solar radiation f to go through water depth. This depth depends on the wavelength of the irradiance and the optical properties of the water. Short wavelengths can go further down, and clear water will let the light go further down.

From the available data, [START_REF] Jerlov | Optical oceanography[END_REF] has constructed the irradiance transmittance of 10 water types and 16 different wavelengths. [START_REF] Paulson | Irradiance Measurements in the Upper Ocean[END_REF] then proposed a method that uses only 'shortwave' and 'longwave' irradiances. They processed the data of Jerlov and provided a set of data with "shortwave attenuation length", "longwave attenuation length", "percentage of shortwave", "percentage of longwave".

The SYMPHONIE model uses the method proposed by [START_REF] Paulson | Irradiance Measurements in the Upper Ocean[END_REF] for parameterizing the light attenuation. By default, the shortwave attenuation length follows the type I of seawater classification [START_REF] Jerlov | Optical oceanography[END_REF]) with a unique value of 23 m over the whole domain. However, the optical properties vary over space and time due to several factors (IOGC, 2000;[START_REF] Saulquin | Estimation of the diffuse attenuation coefficient K dPAR using MERIS and application to seabed habitat mapping[END_REF]: biological activity (mainly the Chlorophyll concentration of phytoplankton), inorganic particulate suspended material (e.g. from sediment resuspension or from dust deposition) and coloured dissolved organic matter (CDOM) that can be produced locally or originate from rivers. In the coastal ocean with sediment transport from rivers, the SWAD is lower than in the open ocean. It also varies over time due to the temporal variation of the sediment transported to the ocean. Previous studies (e.g. [START_REF] Murtugudde | Implementation of quadratic upstream interpolation schemes for solute transport into HYDRUS-1D[END_REF] for the open ocean, [START_REF] Maraldi | NEMO on the shelf: assessment of the Iberia-Biscay-Ireland configuration[END_REF] in the Bay of Biscay/North Atlantic) have evidenced the large sensitivity of simulated temperature and stratification to spatially varying SWAD. We propose to investigate the sensitivity to this parameter in our simulation with the aim to improve the representation of the temperature field with respect to the available observations

II.2.3.2 The light attenuation schemes

The assumption of light penetration from the surface to depth is modeled, following [START_REF] Kraus | Atmosphere-ocean interaction[END_REF], as follows:

𝐼 (%) = 𝑄 2 * t 𝑅 D × 𝑒 & % 3 + E DF+ Eq. 2.33
With:

• 𝐼 (3)
Penetrative solar radiation at depth z • 𝑄 ( solar radiation at the surface, assumed to be made of n radiations i

• 𝑅 5 percentage of radiation i, with ∑ 𝑅 D D = 1
• 𝑧 depth, with the axis at the surface, positive downward • 𝐿 5 penetration length for the radiation i

This equation allows a fraction of radiation (infrared) to be absorbed near the surface while another fraction (blue) can penetrate deeper [START_REF] Kraus | Atmosphere-ocean interaction[END_REF]. In some cases, to simplify the model, n = 1 is used [START_REF] Murtugudde | Implementation of quadratic upstream interpolation schemes for solute transport into HYDRUS-1D[END_REF]. It means that all the radiations unite and share the same attenuation length. In some other cases, to specifically characterize the penetration of different radiations, n = 3 or more can be used [START_REF] Sweeney | Impacts of shortwave penetration depth on large-scale ocean circulation and heat transport[END_REF]. Typically, n = 2 is used for many studies [START_REF] Maraldi | NEMO on the shelf: assessment of the Iberia-Biscay-Ireland configuration[END_REF][START_REF] Ling | A multilevel ocean mixed layer model resolving the diurnal cycle: Development and validation: MULTILEVEL OCEAN MIXED LAYER MODEL[END_REF][START_REF] Zhou | Impact of different solar penetration depths on climate simulations[END_REF].

In SYMPHONIE, n = 2 is chosen. The original equation becomes:

𝐼 (%) = 𝑄 2 * x𝑅 × 𝑒 & % 3+ + (1 -𝑅) × 𝑒 & % 3" y
Eq. 2.34 with:

• 𝑅: percentage of red and near-infrared in the solar radiation (58%)

• 𝐿1: penetration depth of red and near-infrared (set = 0.35m)

• 𝐿2: penetration depth of visible and ultraviolet.

In Eq. 2.31, the first part represents the percentage of energy which comes from the longwave irradiance. Due to low energy, this irradiance can penetrate just a few cm down from the surface and the longwave attenuation length (L1) is set equal to 0.35m for water type I [START_REF] Paulson | Irradiance Measurements in the Upper Ocean[END_REF].

The latter part of the equation corresponds to the irradiance from the short-wave visible light and lower band, like ultraviolet. It is called shortwave attenuation length (L2), or sometimes "Shortwave attenuation depth" (SWAD). SWAD varies over space and time depending on the different optical properties of the water. To deal with it, Jerlov (1968) attempted to characterize L2 in different ocean areas with different water types based on their optical clarity. For example, in water type I, the attenuation length is equal to 23m -this is the parameter value currently used in the Piton et al (2021) configuration. However, since the study area is characterized by strong seasonal activity of river transport near the coast, the use of spatially-varied L2 is expected to produce better results than using a single value.

L2 can be estimated from satellite products. The term "diffuse attenuation coefficient for the downwelling photosynthetically available radiation" (Kd(PAR) or kPAR), which is proportional to L2, is used to mention the ability of light to go through depth.

𝐿2 = 𝑆𝑊𝐴𝐷 = 𝑃𝐴𝑅 = 1 𝑘 678 Eq. 2.35
To calculate kPAR, the kD490 (diffuse attenuation coefficient at 490 nm) is generally used [START_REF] Morel | Examining the consistency of products derived from various ocean color sensors in open ocean (Case 1) waters in the perspective of a multi-sensor approach[END_REF][START_REF] Pierson | Relationship between the attenuation of downwelling irradiance at 490 nm with the attenuation of PAR (400 nm-700 nm) in the Baltic Sea[END_REF][START_REF] Wang | Retrieval of diffuse attenuation coefficient in the Chesapeake Bay and turbid ocean regions for satellite ocean color applications[END_REF][START_REF] Zhou | Impact of different solar penetration depths on climate simulations[END_REF]. kD490 is the attenuation coefficient of diffuse downwelling irradiance at 490 nm [START_REF] Pierson | Relationship between the attenuation of downwelling irradiance at 490 nm with the attenuation of PAR (400 nm-700 nm) in the Baltic Sea[END_REF], a parameter that indicates the turbidity of the water column. In other words, it indicates the ability of the light at 490 nm to penetrate in the water column. Higher kD490 value means smaller attenuation depth, and lower clarity of ocean water. kD490 is acquired by processing data from ocean color satellite sensors such as the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) or the Moderate Resolution Imaging Spectroradiometer (MODIS) [START_REF] Wang | Retrieval of diffuse attenuation coefficient in the Chesapeake Bay and turbid ocean regions for satellite ocean color applications[END_REF]. [START_REF] Morel | Examining the consistency of products derived from various ocean color sensors in open ocean (Case 1) waters in the perspective of a multi-sensor approach[END_REF] proposed an algorithm to calculate the kPAR from kD490 which is applied to optimize the Iberia-Biscay-Ireland configuration described by [START_REF] Maraldi | NEMO on the shelf: assessment of the Iberia-Biscay-Ireland configuration[END_REF]. But this algorithm is valid only to water case I where chlorophyll concentration controls optical properties. For other regions where optical properties are dominated by CDOM, other algorithms are also developed [START_REF] Wang | Retrieval of diffuse attenuation coefficient in the Chesapeake Bay and turbid ocean regions for satellite ocean color applications[END_REF][START_REF] Pierson | Relationship between the attenuation of downwelling irradiance at 490 nm with the attenuation of PAR (400 nm-700 nm) in the Baltic Sea[END_REF].

II.2.3.3 Estimation of attenuation length in Gulf of Tonkin

Saulquin et al. (2013) propose an alternative algorithm to calculate kPAR in the turbidityprevailing region. They suggest 2 different equations for low and high turbidity water:

𝑘 H>< = 4.6051 × 𝐾 𝐷490 6.07 × 𝐾 𝐷490 + 3.2 𝑖𝑓 𝐾 8IJ! ≤ 0.115𝑚 &+ Eq. 2.36a 𝑘 H>< = 0.81 × 𝐾 8IJ! !.K"7L 𝑖𝑓 𝐾 8IJ! > 0.115𝑚 &+ Eq. 2.36b
In this equation set, the first one corresponds to the clear water and the second to the turbid water.

To calculate L2 for the GOT region, we use the monthly means of kD490 for 9 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) from the GLOBCOAST1 project (Jamet and Dessailly, 2016) and distributed by IFREMER. kD490 is calculated using the data from the MEdium Resolution Imaging Spectrometer (MERIS) sensor, installed on Envisat-1 environmental research satellite, using the method described in [START_REF] Jamet | Retrieval of the spectral diffuse attenuation coefficient Kd(Λ) in open and coastal ocean waters using a neural network inversion[END_REF]. Figure II.7 shows an example of kD490 in March and September, 2010. In general, kD490 is higher in March, both at the center of the gulf and in the coastal areas. In September, near the RR, the area with high value of kD490 is advected further offshore. In both months, kD490 is high near Hainan Island, indicating that there is a source of high turbidity water in both seasons. The high turbidity near the Hainan island in Figure II.7 can come from its small rivers or from the vertical mixing created by tides (as suggested for instance by Nguyen et al., 2014). kD490 is interpolated into the simulation grid. Then, kPAR is calculated from kD490 by using Eq. 2.36 and then averaged annually, as shown in Figure II.8. The calculated L2 here is much shorter than in the clear ocean water (14m maximum compared to 23m) with the presence of regional variability. In the coastal area with river activity (bathymetry < 20m), L2 can fall under 4m while in the deepest area it can increase up to 14m. This spatially varying L2 will be used in the reference simulation. The sensitivity of this change will be assessed in Chapter 3.

II.2.4 Carrying out the numerical simulations

The workflow regarding the run and analysis of the simulation is as follows:

• II.2 lists all the multiyear simulations used for the analyses in the next chapters. The simulation contains all the settings as described in II.1, II.2, and with spatially varying L2, as in II.3, called GOT_REF. Three other simulations are performed over the same period:

-GOT_23 is a sensitivity test that we used to assess the impact of the new L2 parameterization.

-a twin simulation without river forcing (GOT_NORIV) to assess the impact of the river runoff on the coastal circulation; all other forcings and parameters being the same as in GOT_REF.

-a twin simulation without tides (GOT_NOTIDE) to assess the impact of tides on the main patterns of the river plume variability; all other forcings and parameters being the same as in GOT_REF.

Besides, two ensemble simulations are performed to assess the impacts of the uncertainty of atmospheric forcing on the ocean variables in general and on the river plume in particular. The details of the ensembles will be given in Chapter 5. To configure these simulations above, several preliminary experiments had been carried out. I spent my first 9 months of my PhD in France, where I ran the model and tested the new L2 on the local cluster at the Observatoire Midi-Pyrénées (NUWA system). The result of this new parameterization will be assessed in Chapter 3. Then, I came back to Vietnam for 1 year. At that time, we decided to work with an updated version of the code with respect to the version of [START_REF] Piton | Du Fleuve Rouge au Golfe du Tonkin: dynamique et transport sédimentaire le long du continuum estuaire-zone côtière[END_REF], which was ~3 years old. However, the connection to France was very slow so I decided to configure and re-run all new simulations on the local cluster at USTH (HILO system) in Vietnam, and redo all the analysis. At that stage, the plume was identified not by using tracers but by using the SSS. When I came back to France for the last year of PhD, we realized that using SSS to calculate the plume area may lead to an overestimation of the plume which will be discussed in Chapter 4. Therefore, we decided to use the passive tracer to get better results. It again took lots of time to rerun all the simulations with two different boundary conditions of the passive tracer. We also ran some experiments on the regional computing center (CALMIP) in Toulouse. In Table II.2 the list may look very concise, but in the end, more than 20 full time tests were performed which consumed more than 650 000 CPU hours, let alone other small tests. It is equal to around 74 years if we only use 1 CPU.

II.3 Observational data sets used to evaluate the simulations

Before analyzing the results, the model is validated with various observational data, including the salinity and temperature data from CTD and satellite, the tide from tidal gauges and satellite altimetry, and the near surface current from the HF radars. In this section, I will present the description of these data. The model validation itself is presented in the next Chapter. We use temperature (T) and salinity (S) profiles measured at 35 CTD stations (Figure II.9B) during the VITEL cruise which took place in July 2014 (Ouillon, 2014). We also use CTD measurements of T, S acquired repeatedly at 10 stations along a 25 km cross-shelf section (Figure II.9C) by Vietnamese and US teams from Center for Oceanography (CFO) and Oregon State University (Rogowski et al., 2019). This dataset, hereinafter referred to as CFO data, consists in 20 timeframes between September 2015 and July 2016.

II.3.1 In situ data

We also compare the model-simulated tidal amplitude and phase with historical tidal measurements from 13 stations (whose location is indicated in Figure II.9D). The dataset stems from the International Hydrographic Organization (https://www.iho.int/).

II.3.2 High-Frequency radar measurements

We use surface velocity data from the high-frequency radar (hereafter HFR) system based on two antennas located at 18.62°N (XUAN site) and 17.47°N (DHOI site) (Figure II.8A) along the coast and operated by the Center for Oceanography, Vietnam Administration of Sea and Islands (CFO, VASI). The data consists of daily maps of zonal and meridional components of the surface current, over the year 2015, built and provided by Tran et al. (2021). As explained in Rogowski et al. (2019), the summer coverage is lesser than in winter because of low sea state conditions; therefore uncertainties on the interpolated velocities are larger in summer. Comparisons with in-situ measurements over a 12-day period indicate a mean bias of 3 cm/s and an RMS difference of 10 cm/s (Tran et al., 2021); we use these values as rough estimates of the data uncertainties. HFR data is representative of currents at 2.4m below the surface (Tran et al., 2021); as a consequence, for the model assessment, we estimate the model current at 2.4m before interpolating it over the HFR grid.

II.3.3 Altimetric data

Tidal constituents computed from satellite altimetric data provide a rich dataset to evaluate the simulated tides offshore. Along-track amplitudes and phases are calculated from the long time series of sea surface height obtained from 20 years of satellite altimetry (TOPEX/Poseidon, Jason 1-2) by using harmonic analysis. We use the dataset described in Lyard et al. (2021). Location of altimetry tracks is shown in Figure II.9D.

The daily L4 product (hereafter AVISO) distributed by the Copernicus Marine Services is used to evaluate the modeled sea surface height (SSH) . It is created with respect to a twenty-year mean from several altimeter missions: Jason-3, Sentinel-3A, HY-2A, Saral/AltiKa, Cryosat-2, Jason-2, Jason-1, T/P, ENVISAT, GFO, ERS1/2. The data can be assessed from the link below: https://resources.marine.copernicus.eu/productdetail/SEALEVEL_GLO_PHY_L4_MY_008_047/INFORMATION For this assessment, the sea level anomaly (SLA) is computed by subtracting the mean SSH over the modeling period (2011-2016) from the daily SSH of the AVISO. Similarly, the mean simulated SSH is subtracted is removed from the daily simulated SSH.

II.3.4 Other satellite data

SMOS (Soil Moisture and Ocean Salinity) L3 v04 provides global sea surface salinity (hereafter SSS) with a resolution of 0.25 degree every 4 days. More information about the SMOS product and data availability can be found here: https://www.catds.fr/Products/Available-products-from-CEC-OS/CEC-Locean-L3-Debiased-v5 OSTIA (Operational SST and Ice Analysis) use data from several satellites to create and provide global sea surface temperature (hereafter SST) in a regular grid at 0.054-degree resolution. The data is available every day. More information about OSTIA can be found here: https://podaac.jpl.nasa.gov/dataset/UKMO-L4HRfnd-GLOB-OSTIA

II.4 K-means clustering analysis

Many different methods exist to analyze the main characteristics of time series, such as statistical, spectral and classification methods. Clustering analysis is a method to discover the structure of data. For example, one has a set of data, each data is described by a set of features. The objective is to group the data into clusters so that the elements within a cluster have a high degree of "natural association" among themselves while the clusters are "relatively distinct" from one another [START_REF] Anderberg | Probability and mathematical statistics: A series of monographs and textbooks[END_REF]. It was used in many areas, from life science (e.g. biology), earth science (e.g. remote sensing), to computer science (e.g. pattern recognition).

There are several algorithms for clustering analysis. For example, [START_REF] Hande | Observed trends in wind speed over the southern ocean[END_REF] use Kmeans algorithm to identify 4 wind regimes over Macquarie Island. [START_REF] Weidle | Validation of strategies using clustering analysis of ecmwf eps for initial perturbations in a limited area model ensemble prediction system[END_REF] applied Complete Linkage algorithm to assess two strategies of generating initial perturbations in a Limited Area Model Ensemble Prediction System. [START_REF] Delaval | A fuzzy classification of the hydrodynamic forcings of the Rhone River plume: An application in case of accidental release of radionuclides[END_REF] use Fuzzy clustering to identify six different hydroclimatological trends from wind and runoff data, then use the results to force a model to examine the variability of the river plume in 6 different conditions. The Self-Organizing Maps method is used to describe the river plume patterns by Vaz et al. (2018) for the Tagus River and Falcieri et al. (2014) for the Po River.

Among different algorithms, K-means clustering analysis (KMA) is a popular unsupervised learning method (Hastie et al., 2001) that allows to classify objects (observations, model outputs) into different groups, given a measure of dissimilarity. The aim of this method is to iteratively identify clusters by their centroid (the means) then minimize the distance between each member of the cluster and the centroid of that cluster. It was first proposed by [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] and is used widely in ocean science. [START_REF] Yi | A clustering analysis of eddies' spatial distribution in the South China Sea[END_REF] used this algorithm to identify spatial distribution of the eddies in the South China Sea. [START_REF] Hänninen | Plastic debris composition and concentration in the arctic ocean, the north sea and the baltic sea[END_REF] applied KMA to define the 4 oceanographically relevant groups characterized by polymer type, size and particle type. Solabarrieta et al. (2015) used it to examine the relationship between wind and surface circulation in the Bay of Biscay. They found that most of the current patterns are related to a specific wind pattern in the study area. Chen et al. (2017) identified the area of Pearl river plume by applying the K-means clustering to a summer climatological turbidity image. Sonnewald et al. (2019) used it to classify dynamically similar regions of the global ocean based on the barotropic vorticity.

In this study, we use KMA to identify the main patterns of the plume and their temporal variations using daily model outputs. To do this, we use the KMA as implemented in the scikitlearn library coded in python (Pedregosa et al. 2011).

The detail process of the KMA is shown in Figure II.10. Firstly, the river plume surface area is identified using a threshold of tracer concentration. Then, it is converted to a binary format, i.e. each point of the model grid is set to 1 (0) for the tracer concentration higher (lower) than this threshold. This is similar to what Falcieri et al. (2014) did for their EOF analysis of the river plume based on SSS. Finally, the KMA is applied to this binary description of the plume. In KMA, the number of clusters is defined by the analyst. Therefore, one problem is that which is the appropriate number of clusters. In some studies, the number of clusters is the authors' choice without further explanation (Chen et al. 2017), or it is chosen after a sensitivity test with the number of clusters (Solabarrieta et al. 2015[START_REF] Hande | Observed trends in wind speed over the southern ocean[END_REF], or with the help of some indicators [START_REF] Hänninen | Plastic debris composition and concentration in the arctic ocean, the north sea and the baltic sea[END_REF], Sonnewald et al. 2019[START_REF] Yi | A clustering analysis of eddies' spatial distribution in the South China Sea[END_REF]. [START_REF] Anderberg | Probability and mathematical statistics: A series of monographs and textbooks[END_REF] and [START_REF] Guanche | Autoregressive logistic regression applied to atmospheric circulation patterns[END_REF] share the same view that the number of clusters must be decided by the users according to their experiences because the results of the clustering analysis is linked with the structure of the data (e.g. the variability of the plume should be explained by the variability of the related forcing conditions). In this study, the result of the KMA, as well as the explanation of choice on the number of clusters, will be presented in Chapter 4.

CHAPTER III: MODEL ASSESSMENTS AND CALIBRATIONS

***

III.1 Introduction

To assess the model's performance, we compared it with various datasets (as described in Chapter 2). These include satellite data (altimetry for tide and sea surface anomaly, OSTIA for sea surface temperature, SMOS for sea surface salinity), in situ data (VITEL, CFO, both with temperature and salinity profiles), and HF radar (for surface current). We also compare our results with the OGCM (the global circulation model) to emphasize the added value of our high resolution regional model compared to a global model. Some of the comparisons have been done by Violaine Piton in her thesis (HF radars, OSTIA, SMOS) with a different version of SYMPHONIE and different light penetration parameterization. Therefore, before doing any analysis, we will assess the model with these data again. Most of the assessments have been published in the article in Chapter 4. These include the assessment of tides with respect to tide gauges and along-track altimetric data, of surface currents by comparing with HF radar data, and of temperature and salinity profiles close to the delta with respect to in situ data (VITEL CTD). Here I only introduced the remaining assessments that have not been published. The pathway shown in OGCM and GOT_REF is more consistent with the cyclonic current in the winter described in previous studies (Ding et al., 2013;Gao et al., 2014;Rogowski et al., 2019;Piton et al., 2021). It is also possible that the SSS is too high in the OGCM which is used for initial and open-boundary conditions, leading to the overestimated SSS with respect to SMOS data in the center of the basin. In summer, the mean SSS values for SMOS, OGCM and GOT_REF are 32.46,32.9,and 32.4, respectively. SMOS shows the low salinity pattern created by the Red River runoff and higher salinity near the southern boundary. This SSS gradient is well represented in GOT_REF, while in the OGCM, the salinity is the lowest in the north of the gulf and the low salinity area due to the Red River and other Vietnamese rivers runoff is not represented. Following a discussion with the MERCATOR-Océan team who produced the global OGCM, we know that the Red River runoff is misrepresented in this configuration, leading to an overestimation of the SSS along the Vietnamese coast.

III.2 Assessment of simulated temperature and salinity in GOT_REF

In GOT_REF, the area associated with the plume shows salinity values that are much lower than those of SMOS. This could indicate an overestimation of the runoff or of the plume offshore and southward extension. However, we believe that this is likely due to increased errors in the SMOS observations in this area instead, since the core of the plume is close to the coast where SMOS data are unavailable. The accuracy of the measurements at the closest point to the coast may not be as good as further offshore. We were surprised by the error patterns provided with the SSS maps (Figure III.1 G, H); these error maps are empirically derived from the observed spatial SSS variability and might be underestimated close to the coast where the SSS data are partially extrapolated (J. Boutin, pers. comm.). In the summer, the spatial variations in OSTIA are smaller. Over the whole area, the temperature varies around 29-30°C, with lower values near the Vietnamese coast south of 18°N and south of Hainan island (Figure III.2D). Compared to OGCM, GOT_REF simulates lower temperatures, which is more consistent with the pattern in OSTIA. In detail, the low temperature area near Hainan island is presented. It seems to be created by tides, as suggested by Nguyen et al. (2014) and Piton et al. (2020). Indeed, this pattern does not appear in the GOT_NOTIDE (not shown). However, the temperature in that area is underestimated by around 2°C, suggesting that the tidal mixing may be overestimated. Overall, performance of GOT_REF is still better than the OGCM, as shown with higher spatial correlation (0. The AVISO data used in this study is computed from multi-year along-track data, and interpolated to the global grid map. Therefore, it can contain errors far from the tracks of the altimeter. Besides, altimetric data are expected to be less accurate in coastal areas because of instrumental errors and larger uncertainties in the wet tropospheric and geophysical corrections (tides, surface waves). In the future, it would be more relevant to use the along-track data. In September, S shows a strong stratification, with values lower than 27 at the surface increasing to over 33 at 5m depth. The model shows the same pattern, with very low bias (0.2 in absolute value) and high spatial correlation (R = 0.91). T is stratified as well. For both T and S, the model shows a slightly stronger stratification, resulting in more saline/colder values in the deepest part of the section. In January, S is mixed over the whole depth with an average salinity of 30 in the first 20km from the coast and is weakly stratified offshore. The model underestimates S in the first 5 km (S ~28) and does not reproduce the stratification offshore. The same behavior is found for T. In April, S is stratified again throughout the section, but S is higher than in September by 1.4 on average. The model shows a similar pattern (R = 0.91), but S is overestimated by 1.2 on average. Similarly, the correlation is high for T (R = 0.96), but the model overestimates it by 2.2°C. Since the strength of the stratification of the GOT_REF and the data is similar and the bias is similar over the whole water column, this overestimation can be related to the uncertainty of the boundary. In July, S decreases significantly throughout the section in both data and model. In the data, the surface layer is about 10m deep, with a mean salinity of ~32 and a minimum value obtained at about 15km from the coast. The model surface salinity is more uniform: the overall correlation is 0.62. In contrast, the model temperature matches well with the data (bias = 0.3°C).

III.2.1.2 Comparison with OSTIA data

Over all 20 stations available, the mean biases (absolute value) between data and model are 0.61 and 0.71°C for S and T respectively, and the mean correlations are 0.81 and 0.72. In April and July, the comparisons suggest that the model runoff is not high enough (as already suggested by the comparison with the VITEL data). In April, the temperature and salinity are higher in the model in the deepest region of the section, which can be due to a bias in the salinity and temperature imposed at the open boundary condition in the south. The vertical mixing seems underestimated also. Overall, despite having biases in some days, the model reproduces accurately the seasonal conditions even in the very shallow area.

III.3 Assessment of the new L2 parameterization

In the initial configuration of Piton et al. (2021), the L2 parameter (shortwave attenuation length) is equal to 23m all over the domain (for more information please refer to Chapter 2, section II.2.3). This is a standard value in the open ocean. However, it should be noted that in the coastal area where the sediment load transported from the river is high, L2 should be shorter than in the open ocean. Therefore, I tried to improve the model accuracy using a new L2 parameterization,which leads to values lower than 4m near the coast to around 14m in the deepest area, as described in Chapter 2. In this section, we will first assess its effect on the temperature and salinity distribution and then check if it can produce better results by comparing them with observational data.

III.3.1 Effect on temperature

We compare the spatially averaged bottom and surface temperature. Indeed, since GOT_23 has a longer L2, it should simulate higher (resp. lower) temperature at the sub-surface and bottom (resp. surface) compared to GOT_REF.

At the bottom, the temperature exhibits a seasonal variation. In the area where the depth is less than 30m, the temperature varies from around 17°C in February to around 30°C in September (Figure III.5A). GOT_23 shows higher temperatures indeed, but the difference between the two runs also has a seasonal variation. From November to February, the difference is very small (< 0.1°C) because the water is fully mixed. From April, the difference increases and reaches the maximum in August (0.8°C on average). In the area deeper than 30m, the seasonal variation of the bottom temperature is smaller, from ~20 to 27°C (Figure III.5B). The difference is again very low from November to February, and it reaches the maximum value in June (0.69°C on average). The seasonality of the impact of L2 on the temperature may be directly due to the seasonal variation of the solar radiation itself; the latter is the largest from April to July (Fig. III9A).

At the surface, the impact of the new L2 value on the temperature is similar between the shallow and deeper regions (not shown) so we just show the comparison of the SST averaged over the whole domain (Fig. III.6). The temperature is lowest in February (~20°C) and highest in August (~30°C). At the surface, most of the time, GOT_23 has a higher temperature than GOT_REF, as expected. From November to January, the difference is less than 0.1°C. Then, the difference increases in March. In April and May, the SST of GOT_23 is 0.3°C lower than the SST of GOT_REF. Overall the impact of the new L2 is larger at the bottom than at the surface. Next, we check the impact of the new L2 on the mixed layer depth (MLD). Here, the MLD is defined as the depth at which the density of seawater is 0.03kgm -3 greater than surface values [START_REF] Zhou | Impact of different solar penetration depths on climate simulations[END_REF]. We computed the mean MLD from December to It should be noted that with a higher value of L2, the mixed layer should be deeper due to the deeper penetration of the irradiance [START_REF] Zhou | Impact of different solar penetration depths on climate simulations[END_REF]. Figure III.7 shows that in both runs the MLD has the same spatial distribution, with slightly higher MLD in the GOT_23. This is also true in winter (not shown) although the difference between GOT_REF and GOT_23 can reach larger values (10m locally). However, there is no direct link between the change of MLD and the change of L2: while L2 in GOT_REF is less than 4m in both the Red River coastal area and near the Hainan strait, the difference of MLD between GOT_REF and GOT_23 is less than 0.2m in the west, while it can reach more than 1m near the Hainan strait (Figure III.7C). This is different from the result of [START_REF] Zhou | Impact of different solar penetration depths on climate simulations[END_REF]. In their study at the global scale, the MLD difference has the same regional patterns as the L2 difference. We propose the following explanation for this different behavior: the impact of L2 on MLD is not significant in the regions influenced by the river freshwater discharge where the MLD is mostly controlled by the low surface salinity. Finally, to check the performance of the new L2, we compared the bias and correlation of GOT_23 and GOT_REF with respect to the VITEL and CFO CTD data. Table III.1 shows that the new L2 improves both the mean bias and the mean of the spatial correlation. Using the new L2 values, the mean bias decreased from 0.82 to 0.61°C for VITEL, and from 0.89 to 0.81°C for CFO. The correlations (R) are improved as well. The comparison with VITEL shows the increase in correlation from 0.57 to 0.71°C, and from 0.57 to 0.74°C for CFO. Figure III.8 shows that GOT_REF indeed produces a better stratification compared to GOT_23 with respect to VITEL profiles. Looking into more details Table III.1, we find that the new L2 decreases the maximum biases for both VITEL and CFO, but it also increases the absolute values of the minimum biases.

Since there is no CTD data available offshore, we can not assess the effect of the new L2 in the deep area. However, the comparison with VITEL and CFO shows that the new L2 produces better results in the shallow area (with bathymetry less than 30m), which is the area where the plume is mostly located.

III.3.2 Effect on heat flux

The total heat flux at the atmospheric-oceanic boundary is made up of four components: shortwave solar radiation, longwave radiation, latent and sensible heat fluxes. The ability of the shortwave radiation to go through water depth is characterized by the shortwave attenuation length (L2). Higher attenuation length allows the irradiance to penetrate deeper and vice versa. In GOT_23, since L2 is higher, the solar component of the heat flux can penetrate deeper, therefore the surface temperature is lower than in GOT_REF. However, since the air-sea heat fluxt depends on the SST, a change in SST also leads to a change in longwave, sensible and latent components. In GOT, the solar flux has seasonal variation, with a minimum value of less than 125 Wm -2 in December and the highest value of more than 250 Wm -2 in May (Figure III.9A). From October to January, the net heat flux is negative for both GOT_REF and GOT_23 (the surface of the ocean is cooling) and the difference of the two simulations is less than 1%. The maximum heat loss occurs in December (-176 Wm -2 for GOT_REF). From February to September, the heat flux is positive for both simulations. While the solar flux reached its maximum in May, the maximum net heat gain occurs in April. In May, the net heat flux decreased from 143 Wm -2 (GOT_23) to 129 Wm -2 (GOT_REF) (10%).

The latent, sensible and longwave heat fluxes from GOT_23 and GOT_REF are shown in Figure III.9B. Among the fluxes, the longwave flux has the smallest variation over time. In GOT_REF, the minimum heat loss due to longwave radiation is 35 Wm -2 (in March) and maximum is 53 Wm -2 (in December). The difference between GOT_23 and GOT_REF is mostly smaller than 1 Wm -2 , except in April and May (~ 2 Wm -2 ). For the sensible heat flux, the maximum heat loss is 42 Wm -2 (in December) while the minimum is in May (6 Wm -2 ). As with the longwave flux, the difference between GOT_23 and GOT_REF is generally smaller than 1 Wm -2 except from April to June. The maximum difference occurs in May (2 Wm -2 ). The latent heat flux contributes the most to the difference between GOT_REF and GOT_23. The maximum difference is 10 Wm -2 , in May. In GOT_REF, the minimum latent heat loss is in April (59 Wm -2 ) and the maximum is in December (202 Wm -2 ).

As a conclusion, we find that the impact of the new L2 parameterization on the heat flux at the air-sea interface is negligible at the basin scale, which is consistent with the weak impact on SST. Maps of the seasonal differences between GOT_23 and GOT_REF do not show significant local differences either (not shown).

III.3.3 Effect on salinity

The RMS of the differences between the daily mean salinity of GOT_REF and GOT_23 over 2011 -2016 are shown in Figure III.10 for different depths. The maximum change of salinity is at the surface, reaching more than 1 near the RR coast, and then decreasing quickly with depth. At 20m, the difference is less than 0.2. There is no significant difference between GOT_REF and GOT_23.

The impact of the new L2 on salinity can be due to several reasons. Firstly, it can result from a change in vertical mixing caused by the change of temperature. Comparisons of density profiles close to the coast (not shown) indeed suggest a change in the strength of stratification (even though the MLD does not differ much) but we still need to do more detailed and quantitative analysis to conclude about this effect. Secondly, it can be linked with changes in evaporation induced by changes on SST. However, the impact of L2 on the latent heat flux (the latent heat flux is proportional to the evaporation in the bulk formulae that we use to compute the air-sea heat fluxes) is very small as discussed in the previous section. Besides, the spatial pattern of the impact on the latent heat flux (not shown) is not consistent with the pattern of the impact on the SSS (Figure III.10). At last, we assume that small perturbations triggered by the change of L2 could be amplified in the coastal area because of flow instabilities (for instance at the front created by the plume), leading to the observed misfit between GOT_REF and GOT_23. This misfit could be seen as the signature of internal variability of the system.

III.4 Conclusion

III.4.1 Model assessment

This chapter presents some results of the model assessment that are not included in [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF] (see Chapter 4).

The model is compared to satellite and in-situ data. SMOS, OSTIA and AVISO products derived from satellite data are used to assess the surface salinity, temperature and large-scale sea level anomaly at seasonal time scales. While the model's surface temperature shows good correlation with the observed data, some misfits between the observed and simulated SSS patterns are found in winter. However, this difference can be due to the error of the SMOS product itself rather than the model error. The quality of SMOS data so close to the coast is questionable here.

The comparison with CFO and VITEL in-situ data shows that the model's skills in representing the T,S vertical gradients as well as their seasonal evolution in the coastal region are reasonably good. In spring, we observe a significant bias in both T and S with respect to the CFO profiles; this might come from errors in the large scale field prescribed at the open boundary condition.. Aside from this bias, nothing indicates a major deficiency in the model or configuration.

However, the comparison with VITEL and CFO salinity profiles suggest a surface freshwater deficit. As suggested in Nguyen-Duy et al. ( 2021) the latter could result from an underestimation of the Red River runoff or an overestimation of the mixing at the river mouths (or in the estuaries) due for instance to uncertainties in tides.

The comparison with HF radar currents, which is not presented in this chapter but in the article of [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF] included in Chapter 4, shows that the monthly variability of the surface currents is consistent with the observations in terms of spatial patterns. In terms of velocity amplitude, significant misfits reaching locally 20 cm/s are observed. They could be due to model errors because of uncertainties in wind forcing or to errors in the HF radar products. In the future, in order to avoid uncertainties due to the HF radar data interpolation onto a regular grid, the model currents should be compared directly to radial velocities.

The evaluation of tides by comparison with altimetric and tide gauges measurements (Chapter 4) for the four main components (O1, K1, M2 and S2) indicates an accurate representation of the tidal signature in sea surface height in the center of the basin (complex error inferior to 3.5 cm). Close to the Red River delta, the comparison with coastal tide gauges leads to larger errors (10.2 cm for K1 and O1) as expected in this area where uncertainties in bathymetry remain.

When comparing the model results with the VITEL CTD profiles, we indeed observed some discrepancies between the deepest CTD depth and the one in the model neighboring points (as illustrated in Fig. III.8). Such discrepancies led us to systematically compare the VITEL profiles with several neighboring points of the model field to check the robustness of our comparisons (not shown).

In this study, we used altimetric along-track tidal constituents and seasonal maps of sea level anomaly. In future studies, sea level anomalies along the track should be explored to estimate the geostrophic component of the coastal current. Depending on the quality of the data so close to the continent, they could also be used to assess the sea level variations in the estuaries and at the river mouth and to analyze the tides, river flow and plume interactions there.

III.4.2 Impact of the new parameterization of the light penetration

In this SYMPHONIE configuration, the spatially varying light penetration depth is calculated from satellite ocean color data to better represent the area under high sediment transport conditions. The model results show that with lower attenuation length, the surface temperature will increase up to 0.5°C locally, while the sub-surface and bottom temperature will decrease (1°C maximum). With lower L2, the MLD is slightly decreased also. The maximum change of MLD happens near the Hainan Strait where there is no river discharge. Due to the increase of the surface temperature, the total heat fluxes can be decreased by more than 10% in spring; the impact of L2 on the heat fluxes is not significant in fall and winter. The heat loss is mostly due to the latent heat flux and the maximum loss occurs in May. Compared to the CTD data, using the spatially varying L2 can create better results than using a single value.

In the framework of this study, we do not use a biogeochemical model. Since the change in temperature and stratification can affect the biological activity, we expect that a change in the L2 can lead to a change in the primary production and the improvement of L2 parameterization could be more meaningful. Also, L2 can vary over time. The calculation with KD490 shows that, in offshore areas, L2 is lowest in the early winter (November, December, January) and highest in the summer. However, in the coastal area, L2 near HaLong bay (near the location of VITEL CTD data) shows the opposite pattern with the highest value in winter and lowest value in summer. Since the difference is large (more than 30%), we expect it can make a significant difference locally in the simulation with high temporal frequency. Therefore, future work could be included to adapt the model to a L2 that varies both spatially and temporally and to investigate its impacts on biogeochemical aspects.

Overall, despite small biases, the model reproduces rather satisfying the conditions for all of the variables that have been checked, and we consider that the model results can be used for further analyses.

CHAPTER IV: VARIABILITY OF RIVER PLUME IN THE GULF OF TONKIN

***

IV.1 Introduction

The Gulf of Tonkin is a shelf sea located east of Vietnam and south of China. In spite of its importance, to date, studies focusing on the plume in the Gulf of Tonkin are still scarce. Therefore, there is a need to better understand the variability of the Red River plume at different scales. In chapters 2 and 3, the method and the validation of the model have been examined. In this chapter, we will attempt to examine the plume variability in the mid-field and far-field regions.

In section IV.2, an article [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF] published as part of this thesis is reproduced. In this paper, firstly, we introduce a short summary of the study area and model configuration, as written in Chapter 1 and Chapter 2. Then, the model is validated with other datasets that have not been introduced in Chapter 3. Finally, the variability of the Red River plume is presented using the model results and K-means clustering analysis.

During the review process of the paper above, the reviewers raised some questions. Not all of the answers led to us modifying the paper but they are interesting concerns and worth looking at. Also, some sensitivity experiments are performed after the paper is published. These discussions will be presented in section IV.3.

IV.2 Paper published in Frontiers in Marine Science

IV.2.1 Summary of the paper

In the paper, the study area and model configuration is described as a short version of what is written in Chapter 1 and Chapter 2. Then, the model is validated with several assessments. Firstly, the first baroclinic Rossby radius of deformation (Rd) is calculated. This is the first time that the Rd is calculated for our area. Then, we assess the general circulation from our model with results from other studies. The mean surface current in the winter (from November to March) is dominated by the gulf-scale cyclonic gyre. From April to August, there is a weak cyclonic gyre in the center of the gulf and an anticyclonic gyre in the south of the gulf. In September and October, the coastal current is southward again from 21°N to the southern boundary where it exits the gulf. This result is comparable with the past studies from several authors.

The next assessments are the comparisons with the satellite and in-situ data. The comparisons with tidal gauges and altimetry data show that the model result performs well in general in representing the sea level amplitude and phase for the main four components (O1, K1, M2 and S2), though still have some errors when compared to the tidal gauges in the delta. These errors can be explained by the uncertainty of the model bathymetry and is comparable with other studies in this area. Compared with HF radar data, the model shows good results, though still has some misfits.

The next section aimed at explaining the variability of the Red River plume. The river runoff is dyed with the passive tracer to better identify the river plume. From several experiments, the plume is defined as the area where the tracer concentration is equal or higher than 7 arbitrary units/m 3 . Then, this criteria is used to calculate the surface plume area. The results show that the plume area follows the same variations as the river runoff, albeit with a time lag. In summer, when the runoff is high, the peak of the plume area occasionally reaches more than 40 000 km 2 (about 27% of the GOT area), while in winter, it sometimes falls below 2 000 km 2 . Due to the interannual variability of the runoff and atmospheric conditions, the plume area varies significantly between different years.

To identify the main spatial patterns of variability, K-means clustering method is applied to daily scenes of the plume that were calculated as described above. The cluster analysis identifies the plume regimes and their periods of occurrence without having to pre-define these periods as one would do when computing seasonal averages for instance. We also apply the classification method to the run without tides to evaluate the impact of tides.

From sensitivity experiments, we have chosen to classify the plume patterns into 4 clusters. The first cluster occurs from November to March. In this cluster, the plume is narrow and contained within the inner shelf (bottom depth < 20m), due to low runoff and downwelling wind (winter monsoon) with an intense downcoast current. The second pattern is mainly observed in the transition period of the monsoon (April, May, October) and is wider than the first pattern; both the relaxation of the winter monsoon and the weaker coastal current allow the plume to spread further offshore than in cluster 1. The third pattern occurs when the summer monsoon is the strongest and the runoff increases. The wind is upwelling favorable; the plumes are advected northward and detach locally from the coast. It is the only pattern where the plumes from the different rivers are disconnected. The final pattern coincides with the highest runoff period and the summer monsoon relaxation. The plumes are connected again in a pool of low salinity waters with the largest coverage, spreading both offshore and southward.

The vertical variability of the plume is examined as well. In winter, the plume is usually mixed over the whole water depth. When the summer monsoon arrives at the end of March, the plume starts to detach from the bottom and spreads offshore as a surface buoyant layer. The plume can detach from the coast near 19°N when a coastal upwelling develops. It may deepen offshore also. In particular, the plume is strongly deepened in the middle of August when trapped in a recurrent eddy near 19°N.

Finally, we tried to identify the "receiving basin" from different rivers. In this study, we define the "receiving basin" for different rivers from the passive tracer content with the same threshold as defined in section 4 (concentration > 7 unit/m 3 ), but this time we apply it to each river system.

The results show that in the low runoff period, the receiving basins from different rivers are very narrow and elongated along the coast, while in the high runoff period, receiving basins of Red River and Southern River are connected.

IV.2.2 Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis

We study the daily to interannual variability of the Red River plume in the Gulf of Tonkin from numerical simulations at high resolution over 6 years (2011)(2012)(2013)(2014)(2015)(2016). Compared with observational data, the model results show good performance. To identify the plume, passive tracers are used in order to (1) help distinguish the freshwater coming from different continental sources, including the Red River branches, and (2) avoid the low salinity effect due to precipitation. We first consider the buoyant plume formed by the Red River waters and three other nearby rivers along the Vietnamese coast. We show that the temporal evolution of the surface coverage of the plume is correlated with the runoff (within a lag), but that the runoff only cannot explain the variability of the river plume; other processes, such as winds and tides, are involved. Using a K-means unsupervised machine learning algorithm, the main patterns of the plume and their evolution in time are analyzed and linked to different environmental conditions. In winter, the plume is narrow and sticks along the coast most of the time due to the downcoast current and northeasterly wind. In early summer, the southwesterly monsoon wind makes the plume flow offshore. The plume reaches its highest coverage in September after the peak of runoff. Vertically, the plume thickness also shows seasonal variations.

In winter, the plume is narrow and mixed over the whole water depth, while in summer, the plume can be detached both from the bottom and the coast. The plume can deepen offshore in summer, due to strong wind (in May, June) or specifically to a recurrent eddy occurring near 19 • N (in August). This first analysis of the variability of the Red River plume can be used to provide a general picture of the transport of materials from the river to the ocean, for example in case of anthropogenic chemical substances leaked to the river. For this purpose, we provide maps of the receiving basins for the different river systems in the Gulf of Tonkin.

INTRODUCTION

River plume can be defined in a general way as the region of the coastal ocean where its properties and dynamics are affected by the river runoff (Horner-Devine et al., 2015). Though the river runoff is small compared to the whole ocean water volume, it can impact both the physics and biogeochemistry of the coastal ocean depending on the discharge, the properties of the ocean area (bathymetry, bottom roughness) and external forcing (airsea fluxes, open ocean influence). Furthermore, rivers carry sediments and anthropogenic contaminants from agriculture and industrial activities. Therefore, there is a need to better understand the fate of the river water from the estuaries to the ocean. It is the first step toward the study of the dispersion of the possible contaminations and toward the design of strategies for monitoring and managing the water quality and the health of ecosystems. These are particularly crucial issues in densely populated areas, such as many deltaic regions of Southeast Asia, including the Gulf of Tonkin.

The Gulf of Tonkin (16.9 • N -21.9 • N, 105.6 • E -110 • E) is a small shelf sea located east of Vietnam and south of China with depth less than 100 m (Figure 1A). It is a meso-tidal region, dominated by diurnal constituents (K1, O1) as described for instance by Nguyen et al. (2014) and Piton et al. (2020). The ocean circulation of the Gulf of Tonkin (hereafter GOT) has been examined by several authors (e.g., Ding et al., 2013;Gao et al., 2014;Rogowski et al., 2019; for a recent review see Piton et al., 2021). Using model results and observational data, they all agree that the winter circulation at the gulf scale is cyclonic and driven by wind. However, in summer, the circulation is not so well explained. Wu et al. (2008) showed that the flow through Hainan Strait has an impact on the circulation in summer: if it is inflow, the circulation is mainly cyclonic and vice versa. From high-resolution model pluri-annual simulations, Piton et al. (2021) found a basin-scale anticyclonic circulation in summer. They also suggest that the surface circulation is mainly ageostrophic, as a consequence from the monsoon wind forcing, except along the Vietnamese coast where the southward coastal current has a dominant densitydriven component.

One of the expected drivers of the dynamics is the freshwater input from the continent. In the GOT, several rivers feed the gulf. The main one is the Red River (hereafter RR) system. It is formed by 3 tributaries that connect at Son Tay (Figure 1B) and then split again into several distributaries. On average, the RR's runoff accounts for more than 60% of the total runoff in GOT (Figure 1C).

In spite of its importance, to date, studies focusing on the RR plume are still scarce. Gao et al. (2013) showed that the coastal plume was found near the northern and western coasts of the gulf in winter while spreading eastward and offshore in summer. Rogowski et al. (2019) analyzed the mean seasonal circulation and suggested that the southwesterly wind direction which is prominent during the summer monsoon is the main mechanism preventing downcoast advection of the RR plume. However, that study did not attempt to explain further the temporal variations of the RR plume. In this study, we propose to take the analysis a step further and examine the plume variability in the midfield and far-field regions (as defined by Horner-Devine et al., 2015) in more detail, using high-resolution simulations and an unsupervised learning method.

The objectives of this study are (1) to propose a method to identify the RR plume in the GOT from the model outputs, (2) to describe its development and characterize its variability at different scales and (3) to attempt to describe the physical processes at work. To do this, we use numerical simulations combined with cluster analysis. In section "Methods, Model, and Data, " the model configurations, data sources and the unsupervised learning algorithm used to classify the main pattern of the plume (K-means) are described. The model is then evaluated against several observational data sets in section "General Circulation and Model Assessment." In section "River Plume in GOT: Identification and Variations of Area, " several methods to identify the plume in the GOT are compared which allows us to select the most appropriate one given our purposes. On that basis, the river plume variability at different time scales is examined in section "Variability of the RR Plume, " illustrating the effect of key physical processes. A discussion on the plume classification and a description of the receiving basins from the different river systems in the GOT are presented in section "Discussion, " followed by conclusions in section "Conclusion."

METHODS, MODEL, AND DATA SYMPHONIE Model and General Configuration

SYMPHONIE is a numerical model that solves the primitive, Boussinesq, hydrostatic equations of the ocean circulation (Marsaleix et al., 2006(Marsaleix et al., , 2008) ) on a curvilinear bipolar (Bentsen et al., 1999). Arakawa C-grid with regular sigma vertical levels. The QUICKEST scheme is used for advection and diffusion of tracers [START_REF] Neumann | Implementation of quadratic upstream interpolation schemes for solute transport into HYDRUS-1D[END_REF], while horizontal advection and diffusion of momentum are respectively computed with a 4th order centered and a bi-harmonic scheme, and vertical advection of momentum by a 2nd order centered scheme (Damien et al., 2017). The k-epsilon turbulence closure scheme is implemented as in Michaud et al. (2012).

Our reference configuration, hereafter called GOT_REF, is an update of the configuration of Piton et al. (2021). The model domain covers the GOT area. Due to the variable horizontal grid, the coastal area near the RR mouths has a fine horizontal resolution of 300 m, while near the open boundary, the grid size can increase to 4500 m (Figure 2A). As the coastal area adjacent to the RR delta is characterized by a complex topography, with many islands and islets, a considerable effort has been devoted to the construction of the model bathymetry and shorelines and is described in Piton et al. (2020). In particular, the bathymetry is reconstructed from GEBCO 2014 combined with other sources and field surveys. The bottom drag coefficient follows a logarithmic law (Blumberg and Mellor, 1987), depending upon the bottom roughness length which is set to 1 mm. The parametrization for the solar penetration depth (as described for instance by At the surface, boundary conditions are provided by an atmospheric model and fluxes of momentum, heat and freshwater are computed using the bulk formulae of Large and Yeager (2004). Operational ECMWF analyses (with a spatial resolution of 1/8 • ) are used to provide 3-h wind, precipitation, solar energy, atmospheric temperature, dew-point temperature, surface pressure.

As the present study focuses on the fate of continental water into the coastal ocean, a specific effort has been deployed on the river runoffs implementation; section "River Configurations" is dedicated to its description. The model is run from 2010 to 2016, starting from the ocean state condition as provided by the global OGCM on 01/01/2010. The time step is set at 2 min. Further analyses are calculated during 2011-2016 (i.e., following a 1-year spinup). The model outputs include daily averaged variables as well as instantaneous fields every 12 h. Unless otherwise specified, both components of the current are detided based on an online harmonic analysis. A summary of the general configuration is available in Table 1.

Two other simulations are performed over the same period:

-a twin simulation without river forcing (GOT_NORIV) to assess the impact of the river runoff on the coastal circulation; all other forcings and parameters are the same as in GOT_REF. -a twin simulation without tides (GOT_NOTIDE) to assess the impact of tides on the main patterns of the river plume variability; all other forcings and parameters are the same as in GOT_REF.

River Configurations

In most modeling studies to date, the RR was set up with only one mouth (or input grid point) using monthly or annual climatological runoff because more realistic runoff at the hydrological station was unavailable (Ding et al., 2013;Gao et al., 2013). In this study and in Piton et al. (2021), the river condition is configured as realistically as possible.

Firstly, the delta is represented taking into account seven input grid points representing the mouths of the main RR delta distributaries (Bach Dang, Cam, Van Uc, Thai Binh, Ba Lat, Ninh Co., Day) (Figure 1B). Secondly, the daily RR runoff is obtained from the National Hydro-Meteorological Service Other Vietnamese rivers (Ma, Yen, Lam) at the south of the Red River delta (hereafter referred to as the 'southern rivers') are also taken into account (Figure 1B). As daily runoffs were not available to us, we prescribe monthly climatological runoffs from NHMS. At the north of the Red River delta, 6 rivers (hereafter referred to as the 'northern rivers') are accounted for (Figure 1B), based on the data given by Gao et al. (2013). In general, the runoff of the Red River system alone accounts for 60% of the total runoff to the gulf, while adding other rivers runoff in Vietnam accounts for around 90% of the total runoff (Figure 1C). In detail, the average discharge for the Red River system in low (December, January, February) and high (July, August, September) discharge period equals to 1632 m 3 /s and 4959 m 3 /s, respectively. For the southern rivers, this value is 365 and 2043 m 3 /s. For the northern rivers, it is 164 m 3 /s and 1103 m 3 /s. It is clearly shown that for the Red River, the ratio between high and low discharge seasons is only 3 times, this ratio is 5.6 times for the southern rivers and 6.7 times for the northern rivers. The lower ratio of the Red River can be due to the presence of several hydrological dams upstream.

In order to simulate the pathways of the river water into the GOT, passive tracers are used. Tracers act as dyes, i.e., they do not affect the dynamics. In total, there are 3 tracers (or three colors) meant to distinguish the inputs from the different river systems. The first tracer is added to the runoff of all the Red River distributaries at constant concentration (100 arbitrary unit/m 3 ). The second one is added at the other rivers in the south and the third one at the rivers in the north of Red River, with the same concentration (100 arbitrary unit/m 3 ). Since tracers are injected with the runoff, they are also submitted to the 1 year of spin-up.

Observational Data Sets Used to Evaluate the Simulation

In situ Data

We use temperature (T) and salinity (S) profiles measured at 35 CTD stations (Figure 2B) during the VITEL cruise which took place in July 2014 (Ouillon, 2014). We further use CTD measurements of T, S acquired repeatedly at 10 stations along a 25 km cross-shelf section by Vietnamese and US teams from Center for Oceanography (CFO) and Oregon State University (Rogowski et al., 2019). This dataset, hereinafter referred to as CFO data, consists in 20 timeframes between September 2015 and July 2016.

We compare the simulated tidal amplitude and phase with historical tidal measurements from 13 stations (whose location is indicated in Figure 3). The dataset stems from the International Hydrographic Organization 1 .

High-Frequency Radar Measurements

We use surface velocity data from the high-frequency radar (hereafter HFR) system based on two antennas located at 18.62 • N (XUAN site) and 17.47 • N (DHOI site) (Figure 2A) along the coast and operated by the Center for Oceanography, Vietnam Administration of Sea and Islands (CFO, VASI). The data consists of daily maps of zonal and meridional components of the surface 1 https://www.iho.int/ current, over the year 2015, built by Tran et al. (2021). The radial velocity measurements are gap-filled, interpolated onto a 6-km rectangular grid and detided as described by Tran et al. (2021). As explained in Rogowski et al. (2019), the summer coverage is lesser than in winter because of low sea state conditions; therefore uncertainties on the interpolated velocities are larger in summer. Comparisons with in situ measurement over a 12-day period indicate a mean bias of 3 cm/s and an RMS difference of 10 cm/s (Tran et al., 2021); we use these values as rough estimates of the data uncertainties. HFR data is representative of currents at 2.4 m below the surface (Tran et al., 2021); as a consequence, for the model assessment, we estimate the model current at 2.4 m before interpolating it over the HFR grid.

Altimetric Data

Tidal constituents computed from satellite altimetric data provide a rich dataset to evaluate the simulated tides offshore. Along-track amplitudes and phases are calculated from the long time series of sea surface height obtained from satellite altimetry (TOPEX/Poseidon, Jason 1-2) by using harmonic analysis. We use the dataset described in Lyard et al. (2021).

K-Means Unsupervised Learning Algorithm for Time Series Pattern Analysis

Many different methods exist to analyze the main characteristics of time series, such as statistical, spectral and classification methods. In terms of classification, there are several methods. For example, Self-Organizing Maps method is used to describe the river plume patterns by Vaz et al. (2018) for the Tagus River and Falcieri et al. (2014) for the Po River. Here, we apply another method, K-means, to analyze not only the plume patterns but the associated forcing conditions.

K-means clustering analysis (KMA) is a popular unsupervised learning method (Hastie et al., 2001) that allows to classify objects (observations, model outputs) into different groups, given a measure of dissimilarity.

The aim of this method is to iteratively identify clusters by their centroid (the means) then minimize the distance between each member of the cluster and the centroid of that cluster. The procedure of this method is made of six steps.

Step 1: Choose the number of clusters to compute.

Step 2: Allocate random numbers as the centroids of the clusters.

Step 3: Calculate the distance from each member to the centroids, using the Euclidian distance (see APPENDIX C for more details).

Step 4: compare these distances, then assign the member to the cluster corresponding to the minimum distance.

Step 5: for each cluster, re-calculate the centroid based on the mean of all members which belong to that cluster.

Step 6: repeat step 3 to 5 until the clusters no longer change.

This method has been used in several past studies in coastal oceanography. Solabarrieta et al. (2015) In this study, we use KMA to identify the main patterns of the plume and their temporal variations using daily model outputs.

To do this, we use the KMA as implemented in the scikit-learn library coded in python (Pedregosa et al., 2011).

In the KMA method, the number of clusters is not fixed automatically, it is chosen depending on the application. The number of clusters should lead to an easy interpretation of the classification, since the objective is to reduce the dimension of the ensemble of scenes to analyze, without eliminating too much variability. In this study, we choose 4 clusters. The method to identify the number of clusters is described in APPENDIX C.

We applied KMA to the "masked river plume." Firstly, plumes are identified using tracers (as described in section "River Plume in GOT: Identification and Variations of Area"). Then, for a given day (a scene), each grid point where the plume is considered to be present is masked with 1, while other points are masked with 0. The results of the classification are shown and analyzed in section "Variability of the RR Plume."

GENERAL CIRCULATION AND MODEL ASSESSMENT

The main objectives of this section are to provide a general description of the patterns of the seasonal circulation at the scale of the GOT from GOT_REF, to assess the consistency of model results with previous studies (or evidence their specificities) and to provide a qualitative assessment by comparing to observations. As a description of the general circulation and model verification (in particular using satellite data) have been provided recently by Piton et al. (2021) from a very similar configuration, we do not provide a detailed analysis. We introduce the section by a discussion of the first Rossby radius, continue with the basin-scale circulation, evaluation of tidal elevation and with the comparison to in situ temperature and salinity data and surface current velocity from high-frequency radars in the area of the RR plume.

First Baroclinic Rossby Radius of Deformation

The first baroclinic Rossby radius of deformation (Rd) sets the scale of mesoscale baroclinic instabilities and, as such, should be resolved by the numerical circulation models (see for instance Greenberg et al., 2007). In river plume dynamics, Rd in the near-field or far-field (coastal current) is used to estimate the relative impact of Coriolis force on the plume dynamics with respect to other forcing (runoff, ambient current).

For instance, Yankovsky and Chapman (1997) showed that for surface-advected plumes, the offshore extension is more than four Rossby radii. In this study, Rd is computed from the mean seasonal field of temperature and salinity from our simulation over 2011-2016. The calculation is made by solving a Sturm-Liouville problem as described for instance in Chelton et al. (1998). We use the method developed by F. Lyard at LEGOS for free surface vertical modes of the internal pressure and described for instance in Nugroho (2017). At depths larger than 60 m, the Rossby radius decreases from 10 to 13 km in April-September (Figure 4A) to 4-7 km in October-March (Figure 4B). In shallower areas, it decreases from 3-8 to 1-5 km. Such a seasonal variation is expected as the stratification is stronger in April-September with higher air temperature. A large spatial variability is also observed; by construction it is correlated to the bathymetric variations. In particular, a striking feature is the eastward extension of the shelf, between the Bach Long Vi island and the coast; the Rd does not exceed 4 km there. Close to the coast in both summer and winter monsoon periods, it can drop to less than 1 km. Moreover, the Rd spatial variability also reflects some patterns of the circulation. For instance, a strong mixing occurs west of Hainan, mostly due to tides, as suggested in the literature (e.g., Nguyen et al., 2014;Piton et al., 2020) and in agreement with the local maximum of the K1 and O1 current in our simulation (Supplementary Figure 1). As a consequence, the Rd is only a few hundred meters west of Hainan Island. On the other hand, the RR runoff leads to an increase of stratification along the Vietnamese coast south of 21 • N which has a clear signature in a locally larger Rd (∼4-5 km) in the October-March period.

Assuming that the effective resolution is roughly between six and ten times its grid resolution, we can evaluate the model's ability to resolve Rd. Rd is well resolved in the northwest quarter of our domain where the mesh is well refined in summer, and along the Vietnamese coast between 19 and 21 • N in winter. However, if we consider that the scales of interest are about times the Rossby radius, then the model resolution is sufficient between 18 and 21 • N.

General Circulation

Mean Surface Circulation

The current in the GOT has been studied with mainly a focus on the northern region, due to the lack of observational data in the south (Wu et al., 2008;Ding et al., 2013;Gao et al., 2014). Recently, thanks to the deployment of a high-frequency radar system, the surface circulation in the southwestern area has been documented (Rogowski et al., 2019;Tran et al., 2021). However, the variability of the general surface and subsurface circulation remains relatively undocumented over the whole gulf and more specifically along the Vietnamese coasts.

Figure 5A shows the mean surface current from November to March from GOT_REF for the 2011-2016 period. The circulation in this winter monsoon period is dominated by the gulf-scale cyclonic gyre, as reported by several authors (Ding et al., 2013;Gao et al., 2014;Rogowski et al., 2019;Piton et al., 2021). The coastal current along Vietnam originates from the merging at 20 • N of the westward flow from Hainan Strait and a coastal current from the northernmost shelf. It flows downcoast down to the southern boundary of the domain. It is the dominant feature of the surface circulation with a mean amplitude of about 30 cm/s. South of Hainan Island, an inflow from the South China Sea is simulated in the deepest part of the area with a decreasing intensity from November to March. This inflow is deflected westward and joins the coastal current south of 19 • N, creating a small-scale cyclonic gyre as in Ding et al. (2013).

In April, the summer monsoon sets in and lasts until August. The southward coastal current weakens and becomes intermittent (Figure 5B). When present, it is deflected eastward at ∼18 • N; it then splits into 2 branches. The first one forms a cyclonic circulation in the central basin. The second branch feeds an intense current south of Hainan that flows out from the gulf. This corresponds to the circulation scheme described by Gao et al. (2014) and Rogowski et al. (2019), but, in our simulation, the current variability is large at daily and interannual time scales.

In September and October (Figure 5C), the coastal current is southward again from 21 • N to the southern boundary where it exits the gulf. It reaches its maximum amplitude (∼35 cm/s) and width. Besides, the circulation is characterized by two inflows with large velocity (∼30 cm/s): one from Hainan Strait and the other from the southern boundary in the deep region.

Horizontal Transport

The horizontal transport displays a similar seasonal cycle as the surface current. From September to March, the general circulation is cyclonic, with a large downcoast transport over the shelf in the west. From April to August, the transport inside the gulf weakens significantly. The main specific patterns are a cyclonic gyre in the central area.

The flow direction during the summer monsoon through Hainan Strait is discussed in several papers (Shi et al., 2002;Gao et al., 2013;Chen et al., 2017). In our simulation, the flux is inflow from September to April and an outflow during the summer monsoon. It varies between -0.18 Sv and +0.15 Sv at 110 E, with significant daily and interannual variability (not shown). Our results however are consistent with the situation described by Wu et al. (2008); in case of inflow the circulation is mainly cyclonic, while in case of outflow the overall circulation is anticyclonic.

Besides, the gulf is fed all year long by an inflow from the southern boundary, which is maximum in summer monsoon period. From September to December, it enters through the east and either penetrates northward in the gulf or is deflected westward and joins the coastal current. The rest of the year it enters through the center, is deflected eastward and generates a large outflow south of Hainan Island.

Model Assessment

Evaluation of Tides

The simulated amplitudes in terms of sea surface elevation for the four main tidal constituents (O1, K1, M2, S2) are shown in Figure 3. O1 and K1 have the same amplitude distribution. The amplitude is largest in the north (100 cm and 95 cm for O1 and K1, respectively) and decreases to its lowest value at the south of the gulf, near the boundary. M2 has two peaks of amplitude. One peak is located in the north (50 cm) and another peak is located at 19 • N (30 cm) along the Vietnamese coast. The lowest amplitude is located near RR mouths. These patterns are similar to the ones found by Nguyen et al. (2014) and Piton et al. (2020). We also examined the tidal currents (see Supplementary Figure 1): the O1 and K1 tidal currents are the strongest west and south of Hainan Island and in the Hainan Strait (more than 50 cm/s). Near the Vietnamese coast, they are weaker and range from to 20 cm/s. They reach a local maximum (∼ 30 cm/s) along the coast at ∼18 • N, 106.5 • E. M2 tidal currents vary between 10 and The simulated tides show remarkable comparisons to alongtrack altimetric data (Figure 3): the root mean squares of the complex errors (i.e., the model-data misfits) over the domain are as low as 2.6 cm, 2.8 cm, 1.3 cm, 0.7 cm for O1, K1, M2 and S2, respectively. These values correspond to 5.0%, 6.0%, 7.4%, and 12.9% of the signal. The model-data misfits are homogeneous in the center of the basin. They are the largest near the coast as expected due to the larger uncertainty of both the altimetry data near the coast and due to the uncertainty of model bathymetry in the very shallow area. These values are comparable to those obtained by Piton et al. (2020) with the T-UGOm tidal model over the same domain (2.6 cm, 3.5 cm, 3.0 cm, and 1.2 cm for O1, K1, M2 and S2, respectively) and to those of Nguyen et al. (2014).

Compared to historical tidal gauges data, the complex error is larger: 18.3, 20.9, 6.2, and 2.3 cm for O1, K1, M2, S2, respectively (Figure 3). There are several reasons for these higher errors. Firstly, the observation period for the tidal gauges data that are available is relatively short (15 days to 1 year) and may not be long enough to be accurate and representative of 'mean' tides. Secondly, the largest errors are observed in the area of the RR delta where the extremely complicated coastline with the presence of thousands of small islands makes the accurate representation of tides a challenge for such a configuration. Compared to Piton et al. (2020), GOT_REF performs better for M2 (6.8 cm) and S2 (3.5 cm), but gives poorer results for K1 (10.2 cm) and O1 (10.2 cm). Future work will be dedicated to tuning the properties of the channels (see section "River Configurations") to better simulate the propagation of the tidal wave in the delta.

Comparisons With High-Frequency Radar Currents

Figure 6 shows the data from the HFRs and from GOT_REF in different seasons in 2015, as well as the maps of correlation between simulated and observed u and v over the whole time series (1 year).

In January, the main observed coastal current is southward, with an amplitude of ∼ 30 cm/s at the coast and decreasing offshore. The model shows the same current direction, but its amplitude is underestimated by 20 cm/s (Figures 6A,B). In April, the coastal current is less intense, but the northward current at 107-107.5 • E is stronger. In GOT_REF, both the southward and northward circulations are observed, with a persistent underestimation of the coastal current regarding the HFR data (Figures 6C,D). In May, the coastal current reverses north of 19 • N in both the observations and simulation. At around 18 • N, the simulated downcoast current is deflected eastward consistently with the HFR observations. However, both at 19 • N and 18 • N the simulated current is stronger than in the HFR data (Figures 6E,F). In July, HFR data show a nearly closed circulation. The coastal current flows southward at 19 • N, then rotates at 18 • N and flows northward and finally joins the coastal current at 19 • N again. In GOT_REF, the circulation shows a more complex pattern, with a northward current (of ∼20 cm/s) at 19 • N and a large temporal variability (Figures 6G,H). In October, HFR data show a similar current pattern as in January, with a strong southward coastal current (∼40 cm/s). The model shows very good consistency, with the strongest velocity near 18 • N (Figures 6I,J).

The temporal correlation between simulated and observed daily fields of u and v is higher for u offshore (∼0.5) (Figure 6K), while the correlation for v is better near the coast (>0.6) (Figure 6L). The larger misfits between the model and data at ∼19 • N seem to come from a larger spatial variability of the simulated u component (not shown) than the observed one. In particular, as discussed in section "Variability of the RR Plume, " this area is characterized by some eddy activity in summer which may not be resolved by the HFR observations or may be out of phase with the observations.

Temperature and Salinity Profiles From in situ Measurements

Figure 7 shows the salinity and temperature profiles from the VITEL campaign and the GOT_REF simulation. Although it is a very shallow area, the bathymetry of the model is quite accurate compared to the depth of the CTD data. The largest bathymetry misfit is reached at stations 26 to 30 where the model is too shallow by ∼10 m. Since the campaign took place in summer (July), there is a strong stratification in both salinity and temperature profiles. At the surface, the observed salinity can be as low as 15, while at the bottom it is around 34. The mean value at all points is 28.7. The model performance is good, with a correlation R = 0.88 computed from all profiles, but the mean value is overestimated by 0.7. The model overestimates the surface salinity at stations located in the shallowest area. This suggests that the river runoff is underestimated or that the simulated mixing in the estuaries is too strong, possibly due to local errors on the tidal representation.

The observed temperature goes up to 33 • C at the surface and 27 • C at the bottom, with a mean value of 30.3 • C. The mean bias of the model is low (0.2 • C), although the temperature is overestimated at both the surface and bottom at some stations. Overall, R is 0.72.

The comparison with the T, S profiles from the CFO campaigns is not described in detail as we get similar results to those described by Piton et al. (2021) in a similar configuration of the model. For all the 20 sections available, the mean biases (absolute value) between data and model are 0.61 and 0.71 • C for S and T respectively, and the mean correlations 0.81 and 0.72. The comparisons suggest that the model runoff is not high enough (as already suggested by the comparison with the VITEL data), or the background salinity provided by the boundary condition is too high. The vertical mixing seems underestimated also. Overall, despite small biases, the model reproduces accurately the seasonal conditions even in the very shallow area.

RIVER PLUME IN GOT: IDENTIFICATION AND VARIATIONS OF AREA

This section will first review some past studies on plume identification then explains our method to identify the plume locally in the GOT. Then, the resulting variations of the plume area in the period 2011-2016 are analyzed.

A Brief and Non-exhaustive Review of Past Studies

There are no clear criteria to identify a river plume in general. Usually, such criteria involve sea surface salinity (hereafter SSS) because the river plume has a lower surface salinity than the surrounding ambient waters. However, due to different environmental conditions, the reference salinity is difficult to set in a general enough manner. In some studies, the authors may choose a suitable value based on their knowledge and experience; however, the choice is seldom justified in a rigorous way. This may lead to a situation where different authors choose different values for the same plume, depending on the observational or modeling approach and the objectives. For example, for the Columbia river plume, Liu et al. (2009) Other authors defined other criteria to identify the plume. Otero et al. (2008) use the highest horizontal gradient of mixed layer depth to classify the river plume area in the Northwest Iberia and then use the SSS value of 35.6 that coincides with this maximum gradient. Another approach to identify river plumes is based on ocean color satellite images that reveal turbid water masses rich in Chlorophyll, suspended matter, or colored dissolved organic matter (e.g., Chen et al. (2017) for the Pearl River).

To our knowledge, there have been very few studies dedicated to the RR plume, so there is no published reference ambient salinity value. Rogowski et al. (2019) locate the RR plume using monthly means of ocean color maps from MODIS. In tropical areas, the presence of clouds is however a limitation to the use of optical products to investigate the sea surface variability at time scales shorter than a month.

Stratification Index and SSS Threshold

To identify the RR plume, different methods are in consideration in this paper. Firstly, we apply a method based on the Stratification Index (hereafter SI), assuming that the river water forms a buoyant layer over the ocean water. This index was first proposed by Hansen and Rattray (1966) for estuary classification. It is defined as the relative difference between the surface and bottom salinity (eq. 1). In this paper, we tested this method with a threshold of 0.1. If SI equals or exceeds 0.1, the water column is defined as belonging to the river plume. The time series of the daily plume area for all the rivers in GOT based on SI is shown in Figure 8A. As expected, the river plume area varies with the river runoff. In July-August (high runoff period), the plume area reaches its highest value, while it is lowest in December-January (low runoff period).

SI =

salinity bottom -salinity surface vertical mean salinity

However, there is one major drawback to the SI method. Vertical salinity sections indicate that the plume in the GOT consists of 2 zones. One is the very shallow coastal area where the plume is mixed over the entire water column; in the second one, farther offshore, the plume is detached from the bottom and stays on top of the high salinity water. SI being based on the difference between the salinity at the top and the bottom of the water column, it cannot detect the plume zone in the shallow area where it is totally mixed. In the low runoff period, as we will show in section "Variability of the RR Plume, " the downcoast current and the northeasterly wind make the plume stick to the coast; in such a case, the whole plume is well-mixed and SI is not helpful since the calculated area is close to zero (Figure 8A).

To deal with that problem, in a second step, we considered using SI as a proxy to find the suitable SSS that can be used as the plume indicator. After several experiments, we selected the value of SSS = 30. In high runoff period, this criterion defines a plume with a similar pattern offshore to the one deduced from the SI FIGURE 8 | (A) Total runoff (in m 3 /s) and plume area (in km 2 ) for all the rivers in the Gulf of Tonkin. The plume area is calculated using the Stratification Index (blue), SSS (red) and tracer concentration (green). (B) Spatially averaged SSS over the whole domain (in blue) from the simulation without river (GOT_NORIV), and difference of river plume area between SSS and tracer methods (in red). (C) Daily runoff and plume area for the Red River and Southern rivers.

criterion. In the low runoff period, this discriminant value also captures the narrow plume area that runs along the coast (Figure 8A).

Using SSS, the envelope of river plumes created by all the rivers in the GOT was generally identified. However, this "stratificationbased SSS" method has in turn two drawbacks. Firstly, it cannot distinguish the plumes created by the different rivers that feed the GOT. Secondly, the SSS criterion may be slightly biased at seasonal time scales by the effect of precipitation. Indeed, the precipitation has seasonal variations due to the monsoon system. The average precipitation from January to March is 1.8 mm/day, while it is 5 times higher (9.4 mm/day) from August to October (calculated from ECMWF daily precipitation over 2011-2016).

Passive Tracers

As an alternative method, we use passive tracers that behave in the model as any buoyancy-free particle or passive chemical from the rivers (as explained in section "River Configurations"). This is particularly relevant to identify the coastal area influenced by the river input of contaminants and to distinguish the freshwater input from other sources than the river such as advection by coastal current and other coastal runoffs [see for instance or to study mixing between the plume and ambient waters [e.g., Vlasenko et al. (2013) for the Columbia River]. As for the SSS criteria, the choice of a tracer concentration threshold to define the plume is author-dependent. In Vlasenko et al. (2013), the river plume is identified as the area where the concentration exceeds 10 units/m 3 while in Wang et al. (2014), the threshold is 5 units/m 3 . In our study, several experiments led to identify the river plume as the area where the concentration exceeds 7 unit/m 3 , which best fits with the area identified by the criterion on SSS in the low runoff period (dry season).

As indicated in section "Observational Data Sets Used to Evaluate the Simulation, " 3 tracers were added, respectively to the RR distributaries, the southern rivers and the northern rivers. Figure 8A shows the plume area for all the rivers in the GOT calculated by summing the three tracers' concentrations. The plume areas computed from the SSS and tracer criteria have a similar trend in both the low and high runoff periods. However, in the high runoff period, especially during the peaks in 2012 and 2013, the plume area calculated from the tracer is much lower than from the SSS. This difference can go up to over 20,000 km 2 (Figure 8B). Figure 8B also shows the SSS of the twin simulation without any river runoff (GOT_NORIV). Even with no river, the SSS shows significant variations between 31 and 33. The low sea surface salinity often happens in summer, which is in the rainy season. It is minimum in 2012 and 2013, which are also the years when the runoff reaches the highest peaks (Figure 8A), hinting at a correlation between local precipitation and runoff in both years. The SSS drop in the GOT_NORIV run has a similar trend than the difference of plume area calculated by the two methods, which suggests that this difference could be indeed due to precipitation.

In conclusion, for our objectives, the method based on the passive tracers seems the most appropriate to help identify the RR plume among the three methods investigated in this study. It can capture the plume in both the highly stratified zone and the shallow mixed zone. Besides, it can distinguish the plumes created by different rivers (as discussed in section "Discussion").

Temporal Variability of the Plume Area

South of 20 • N, the RR plume is quickly joined by the plumes of the southern rivers (Ma, Yen and Lam, Figure 1), creating a unique buoyant plume which extends southward along the coast most of the year. In this section and in the following we therefore analyze the variability of the resulting plume from the RR and the three southern rivers together. Section "Discussion" will discuss the receiving basins by distinguishing the RR from the other rivers, thanks to the multi-tracer approach.

Figure 8C shows the evolution of the cumulated runoff from the Red, Lam, Yen and Ma rivers and the corresponding plume area from 2011 to 2016. As expected, the plume area follows the same variations as the river runoff, albeit with a time lag. In summer, when the runoff is high, the peak of the plume area occasionally reaches more than 40,000 km 2 (about 27% of the GOT area), while in winter, it sometimes falls below 2,000 km 2 . Due to the interannual variability of the runoff, the plume area varies significantly between different years. The river plume area peak in 2013 (41,000 km 2 ) is nearly twice as large as the peak in 2015 (22,000 km 2 ). Both the total discharge and the plume area undergo a strong variability at shorter time scales of a few days as well.

Figure 8C also shows two more characteristics of the plume variability. Firstly, a higher runoff does not ensure a larger river plume area. In 2012, the runoff peak is higher than in 2013, while the plume extension in 2013 is larger. Similarly, although the rainy seasons in 2014 and 2015 are characterized by similar runoffs, the plume area in 2015 is much lower than in 2014. Secondly, there is a time lag between the runoff peak and the peak of the plume area. This lag is also described in Rogowski et al. (2019). Depending on the year, it can be up to 1 month (in 2013). In 2015 and 2016, the plume area is relatively small until July and reaches its peak value rapidly within a few days. In short, both the time lag and plume area do not appear directly correlated to the runoff intensity (i.e., a higher runoff will not necessarily create a larger plume area about 1 month later). Such a variability evidences the fact that river runoff intensity is not the only factor driving the plume size and variability. Wind variability and its impact on the surface circulation and mixing is likely to be another driver, as we show in the next section.

Figure 8C displays as well some interannual variability in both the runoff and the plume area variations. The years 2015-2016 were marked by an intense El Nino event. Generally, El Nino events lead to decreased rainfall, drought and saltwater intrusion in deltas in Vietnam (Sutton et al., 2019). Figure 8C does not evidence a drastic drop in the 2015-2016 runoff with respect to 2014 and 2011. However, the SSS signal in the GOT_NORIV run is significantly different in 2015-2016 with respect to the previous years (Figure 8B): the mean SSS stays above 32 and the seasonal signal is inverted with a maximum SSS in summer. El Nino is likely to influence the coastal circulation, through alteration of wind and air-sea fluxes of freshwater and heat. The specific influence of ENSO on the coastal plume will be investigated in future studies as we focus in this paper on seasonal and shorter time scales.

VARIABILITY OF THE RR PLUME Spatial Patterns of the Plume

In the previous section, the river plume in GOT has been identified using passive tracers. We have documented the time variations of the total surface area and compared them to the runoff variations. Now, the seasonal variability of the plume spatial patterns is examined using KMA.

Figure 9 shows the spatial distribution of frequency of occurrence of the plume within each of the 4 clusters calculated from the plume area from 2011 to 2016, and the temporal evolution of the clusters, revealing both seasonal and interannual variations. The main variables of interest (SSS, wind, surface current, SSH) are averaged over the period corresponding to each cluster and are shown in Figure 10.

The first cluster usually occurs from October to March and appears on 992 days (45.3%) in 6 years. In this cluster, the plume is very narrow and is mostly confined to the shallow area (depth < 20 m). 75% of the time that the plume occurs, at 19 • N, it extends to 105.8E (∼20 km from the coast) while 25% of the plume extends to 105.85E (∼25 km from the coast) (Figure 9A). If the runoff (or other forcing) undergoes strong fluctuations, the plume area should also vary. In this cluster, the difference between the spatial occupancy of 25% and 75% is small: it means that the forcing condition is relatively stable. The simultaneous occurrence of several conditions can explain the relatively small extent of the plume in this cluster. Firstly, it happens in the lowest runoff period (average discharge over the cluster equal to 2197 m 3 /s). Secondly, the wind is strong and from the north-east (winter monsoon, Figure 10A), which is downwelling-favorable. At the GOT scale, the mean wind velocity is 6.4 m/s; it is stronger than in any other cluster, even near the coast (it reaches ∼ 5.7 m/s on average inside the red box shown on Figure 10A). The coastal surface current is southward, with a mean speed larger than 20 cm/s between 20.5 and 17.5 • N (Figure 10E); consistently with the strong monsoon wind, the current also reaches its largest intensity in cluster 1. All these conditions favor the low salinity (SSS < 30) water to be confined to the coast and to extend southward all along the coast (Figure 10I). We also observe a crossshore gradient of SSH with an elevation larger than 5 cm with respect to the basin-averaged SSH (Figure 10M). Here, the coastal current is fed by buoyant river waters which contribute to the geostrophic component of the current. Indeed, in the GOT_NORIV simulation, the surface current in cluster 1 is 10 to 15 cm/s weaker; the SSH is 1 to 4 cm lower (not shown). Similarly, Piton et al. (2021) estimate that the geostrophic contribution to the downcoast surface current reaches up to 60% in December-February.

The second cluster appears on and off throughout the year, with a slightly greater rate of occurrence in April and September (usually before the occurrence of cluster 3 and after the one of cluster 4), and accounts for 518 days (23.6%). It happens mostly in late spring and early winter, i.e., during the seasonal transition of the monsoon. This suggests that cluster 2 represents Nguyen-Duy et al.

Red River Plume Variability a transition regime for the river plume. In this cluster, the river plume extends both further offshore and further southward compared to cluster 1. 75% of the time, the plume extends to 105.85 • E (∼25 km from the coast) while 25% of the time, it extends to 106 • E (∼40 km from the coast). The low salinity strip defined by the 30 isohaline is about twice as wide as in cluster 1 (Figure 10J). 50% of the time, the plume reaches 18 • N, which is the same as in cluster 1, while 25% of the time, it extends as south as 17.1 • N (Figure 9B). The higher spatial coverage of the plume area in this cluster can be explained by the higher runoff (4020 m 3 /s) than in cluster 1 (2197 m 3 /s). Also, in this cluster, the wind direction is from the southeast, therefore different with respect to cluster 1 (Figure 10B). The wind speed is weak, with a mean speed of 4.3 m/s near the coast (4.7 m/s at the gulf scale).

The downwelling effect is relaxed and the light plume water can spread seaward. The coastal current is still flowing southward but is weaker than in cluster 1. As in cluster 1 and cluster 4 (see below), it is locally intensified, around 106.5 • E, 18 • N. The tidal current is intensified there as well (Supplementary Figure 1), probably because this is an area where the shelf is thinning and the coastline draws a cape. When present, the coastal current may be intensified there, through the same processes as for the tidal current. Another assumption is that interactions between tides and the coastal current occur there and result in an intensification of the local circulation. The third cluster, which happens primarily in June and July (315 days, 14.4%), has a different shape with respect to clusters 1 and 2 (Figure 9C). In this cluster, the plume is advected northward and seaward. This is also the cluster in which the river plumes from the various rivers are disjoint. It corresponds to the period when the runoff increases with the mean value equal to 4631 m 3 /s and the summer monsoon wind reaches its strongest intensity. The mean gulf-scale wind speed is 5.7 m/s (5.6 m/s inside the box). It blows northeastward therefore corresponding to upwelling conditions. The surface coastal current is not present anymore south of ∼18.7 • N; at 18 • N the surface current is oriented seaward, with weak speed (∼10-15 cm/s). North of 18.7 • N, the coastal current is reversed, flowing northward, with a weaker speed than in cluster 2 (∼15 cm/s), but locally intensified by the river runoffs. Although the average runoff is larger than for cluster 2, the very low salinity strip (SSS < 25) is smaller than in cluster 2, suggesting that strong mixing is taking place, either vertically due to wind-induced turbulent kinetic energy or horizontally. Figure 9C shows that 5% of the time, the plume reaches the 40 m isobath, and indeed the salinity over this area is decreased with respect to cluster 2, supporting the hypothesis of lateral mixing and/or stirring. The involved processes are assumed to be at daily timescales or shorter (inertial, tidal time scales) and are discussed in section "Variability of the Plume Dynamics at Daily Time Scales."

The fourth cluster exhibits the largest spatial coverage and occurs mostly in August and September (367 days, 16.7%), that is at the peak of the high runoff season (the mean runoff over this cluster reaches 7336 m 3 /s). In this cluster, the plume extends the farthest offshore. At 19 • N, it can extend to more than 100 km from the coast. 75% of the plume extends to ∼55 km from the coast, while 25% of the plume extends to 100 km from the coast. 5% of the plume can extend as far east as 107.0 • N, which is 145 km from the coast (Figure 9D). In this cluster, the mean wind pattern is still upwelling-favorable, i.e., capable of driving the plume offshore (Figure 10D). However, the wind speed is less intense: 4.5 m/s at the gulf scale and 4.0 m/s inside the box only. It is the period with the lowest wind speed and the weakest coastal current; the latest is downcoast again. Its speed does not exceed 10 cm/s (except around 106.5 • E, 18 • N where it is locally accelerated up to more than 20 cm/s as in clusters 1 and 2, Figure 10H). The plume is surface-advected, which is verified by the analysis of the plume thickness in the next section. The large runoff during the August-September period is responsible for the very low salinity (SSS < 25) simulated along the coast between 21 • N and 19.5 • N. As for cluster 3, the salinity over the whole shelf in the western half of the GOT is decreased in this cluster, with respect to the previous one, indicating a larger extent of the river water influence than the one formalized by the plume definition that we adopted in section "River Plume in GOT: Identification and Variations of Area." In other words, even highly diluted in the ocean waters, the river water is still influencing the ocean SSS. The comparison between the reference run and the run without runoff (GOT_NORIV) confirms that the runoff influences the SSS all over the gulf in clusters 2, and 4 (Supplementary Figure 2), with the largest influence in cluster 4 (more than 1 west of 107 • N and north of ∼20.5 • N). The SSH shows a small elevation between the center of the basin and the Vietnamese coast, consistently with a weak downcoast circulation (Figure 11P). In GOT_NORIV, the SSH east of 107 • E is decreased by 1 to 2.5 cm with respect to REF (not shown), while the surface current only differs from the REF run locally and the differences are smaller than for cluster 1 (and to a lesser extent for cluster 2), therefore suggesting a weak contribution of the riverine waters to the circulation in late summer.

We estimate the export of riverine waters by computing the daily fluxes of dye concentration (for the RR and southern rivers) through three sections (Figure 2A). We then average the fluxes over the periods of the cluster to obtain a mean daily flux for each cluster. We found that the highest offshore flux (section 2) occurs in cluster 3 (∼402,000 unit/s), while the largest downcoast flux (section 1) is in cluster 2. Overall the largest exports are across section 1 in clusters 2 (473,532 unit/s) and 4 (431,408 unit/s). The flux across section 3 is much smaller in all the clusters; it is the largest (northward) in cluster 3 (76,000 unit/s).

Plume Thickness

We define the plume thickness as the maximum depth at which water is found with a tracer concentration greater than 7 unit/m 3 . It is computed daily. Figure 11 shows the time evolution of plume thickness and temperature at 5 m along the 19 • N section. For the sake of clarity, we show years 2013 and 2014 only, representative of high and low discharge conditions respectively, but similar conclusions would be drawn for the other years (Supplementary Figure 3).

As we have seen in the previous section, from September to March, the plume is narrow and elongated along the coast due to the winter monsoon downwelling winds. During that period the plume is bottom attached, i.e., it is filling the whole water column with a thickness of the order of 10 m at 20 • N and 15 m at 19 • N (Figure 11A). Several factors can explain such a thickness: first the downwelling winds tend to create a convergent flow onshore; secondly the vertical mixing may be enhanced by the bottom-generated turbulence resulting from the strong coastal current (Wiseman and Garvine, 1995). Wiseman and Garvine (1995) also suggest that onshore Ekman flow may result in dense ocean water overriding buoyant plume water, therefore creating instabilities and mixing. From May on, the monsoon changes to southwesterly wind, driving the plume northward and detaching it from the coast and from the bottom. We observe that events where the plume detaches from the coast at 19 • N coincide with the temperature at the coast colder than seaward (Figure 11B). This supports the assumption that the southerly monsoon winds generate a coastal upwelling. The offshore extension of the low-temperature signal is less than 40 km. This summer upwelling was also simulated by Gao et al. (2013) for 2007, and is likely to be responsible for the detachment of the plume from the coast. When detached, the plume thickness increases from ∼4 m onshore to about 10 m seaward. Sometimes in May-August, the plume is thickening offshore to 12-15 m (Figure 11A). That deepening may be due to strong winds events that increase the mixing offshore in the upwelling region or to the increase of river discharge.

In mid-August, we observe a strong plume deepening. In 2013, the plume is still attached to the coast, contrary to 2014 (Figure 11A). This difference can be due to two reasons. Firstly, the runoff in 2013 is higher than in 2014. Secondly, the wind in 2014 is stronger than in 2013 (not shown). However, both deepening events appear to be linked with a seasonally recurrent anti-cyclonic eddy developing near 19 • N, which is examined in more detail in the next section.

We also look into detail at the temporal evolution of the tracer concentration over the water column at some points near the coast (Figures 11C-E). After a peak discharge in summer, from October until November or December (depending on the year), the tracer is rapidly mixed down to the bottom (∼25 m) at a low concentration (3-5 unit/m 3 ). Then, from January to March, the concentration gets close to zero, except for some sporadic events of a few days.

Variability of the Plume Dynamics at Daily Time Scales

The cluster analysis allowed to identify clear seasonal patterns of the plume variability, but the analysis of the plume thickness variability highlighted a strong variability at daily time scales as well. In particular, the plume is observed to spread seaward in spring and summer, which is synonymous of export of fresh water and riverine materials from the coastal zone to the interior of the GOT. Those sporadic events are represented by the 5% of plume occurrence in Figure 9. Figures 12A-C illustrates some of these events for three dates in June 2012. They show that the river waters spread within the whole basin, with a low concentration though. They extend as far east as Hainan coasts. As the forcings (wind, runoff, tides) are also strongly variable, it is not possible to disentangle their respective impacts.

Besides, the plume seems to be shaped by small-scale circulation patterns over the whole shelf, suggesting that mesoscale or submesoscale activity strongly influences the horizontal transport as well as the dilution.

In particular, at 19 • N, we observe a recurrent anticyclonic eddy as mentioned in the previous section. Figures 12D,E shows the relative surface vorticity and the plume thickness on and 15 August 2014. The eddy is depicted by the minimum of negative surface vorticity. It appears clearly that the plume is deepening at the center of the eddy, with a difference of thickness reaching 8 m between the edge and the center for instance on 15 August 2014. In our simulation, this eddy generally happens in August, when the wind direction is southwest, which is upwelling-favorable, and the runoff is high. Then, this eddy disappears when the coastal southward current develops again and/or the wind is not favorable anymore. Its lifetime varies from a few days to ∼15 days; its diameter is about 50 km, i.e., much larger than the first Rossby radius. There is no mention about such an eddy from the HFR analysis (Rogowski et al., 2019;Tran et al., 2021); the availability and resolution of the HFR measurements may be too limited to depict such a pattern. The formation mechanism of this eddy is unclear. The comparison with GOT_NORIV shows that most of the time, without the river, some vorticity gradients can be depicted, suggesting a weak anticyclonic structure which vanishes on a shorter time scale (for instance in 2013, 6 days in GOT_NORIV instead of 8 days in GOT_REF) (Supplementary Figure 4). The anticyclone seems strongly connected to the Lam river plume. Bottom topography may be another forcing factor. Indeed a small island is present at 106 • N and is likely to influence the river flow. Both the buoyant Lam river plume and the island may generate instabilities of the coastal flow that could lead to the eddy development.

Impact of Tides

The influence of tides in shaping the far field plume is investigated by comparing the reference simulation GOT_REF with the simulation without tides (GOT_NOTIDE). We first compute the total area of the plume from all the Vietnamese rivers (including the RR). On average, the plume in GOT_REF is 5% larger than in GOT_NOTIDE, with the largest differences found in the high discharge period (not shown). Significant differences (∼10%) are observed during the summer of 2011 and 2012, with the plume being larger in the GOT_REF run; over the remaining period, the differences are weaker.

We performed a KMA analysis on the GOT_NOTIDE outputs for 4 clusters; the resulting clusters are very close to the ones from GOT_REF, both in terms of spatial structures and temporal distribution (Supplementary Figure 5). The main differences are found in clusters 3 and 4, that is during the summer monsoon and high runoff period. Without tides, the plume spreads more to the north. In cluster 4, the plume spreads a bit further to the south without tides.

Besides, of equal importance are the tidal effect on vertical mixing (already inside the estuaries) and frontal activities (Guarnieri et al., 2013). The impact of tides on the vertical mixing is investigated by comparing the tracer concentration between GOT_REF and GOT_NOTIDE (Figure 13). The tidal influence differs depending on the area: along the coast, just north of the mouth (close to ∼21 • N), it is larger without tides at both surface and bottom (>4 unit/m 3 ). Everywhere else, the bottom concentration is larger with tides (>0.5 unit/m 3 ). At the surface, the concentration is smaller at the mouths (>4 unit/m 3 ) and larger in the plume (0.5-1 unit/m 3 ) in the run with tides. Over the shelf, the impact is small; in clusters 3 and 4 we observe a larger surface concentration (>0.5-1 unit/m 3 ) in the run without tides at ∼ 107-108 • E and ∼ 20 • N. These differences are small but significant with regard to the concentrations in GOT_REF (Figure 13): the values are 5-50 unit/m 3 (resp. 1-20 unit/m 3 ) in the plume and 0.5-10 unit/m 3 (resp. 0.1-1 unit/m 3 ) over the shelf west of 108 • E at surface and at the bottom, respectively. We interpret the results as follows: tides enhance the vertical mixing, which explains the larger concentration at the bottom in GOT_REF. However, tides also enhance the export of riverine waters offshore at the mouth, hence the slightly larger surface concentration values in GOT_REF in the plume and the smaller concentration close to the mouth. The plume reaches the small coastal area north of the RR mouths only in the absence of tides as mentioned previously, so the surface and bottom concentrations are larger in GOT_NOTIDE.

DISCUSSION

In this section, we summarize our findings, presenting them with paradigms found in the literature on plume dynamics (classification) or on environmental applications (receiving basins).

Classification

Several attempts have been made in the past to classify buoyant coastal runoffs with the objective of deriving general dynamical properties and to compare different plume systems. Garvine (1995) introduced the Kelvin number (K) which is defined as the ratio of the typical cross-shore scale of the buoyant runoff (usually taken as the estuary mouth width) over the local internal Rossby radius. K measures the importance of rotation: for K < 1, the influence of earth rotation on the plume dynamics is not significant and the river outflow forms energetic jets whose direction is controlled by local bathymetry and coastline (Hetland and Hsu, 2013). When K > 1, the Coriolis effect is dominant and the plume flows downcoast forming a coastal current; for large K, the downcoast flow may be insensitive to wind-driven motion in the opposite direction (Wiseman and Garvine, 1995). Very close to the RR mouths, the internal Rossby radius is between 1.5 and 5 km. Given the typical width of the distributaries mouths in our configuration (from 0.6 to 2.5 km), the Kelvin numbers vary between 0.5 and 1. Therefore in theory, bulges could be formed and the impact of the Coriolis force should be small. Close to the mouth, the currents for each cluster do not show evidence of bulge formation. The outflow from the Cam river is very much constrained by the bathymetry and coastline. At the Van Uc and Ba Lat mouths, the outflow forms a kind of jet. In clusters 1 and 2, the plume is deflected downcoast and merges into the strong (wind-generated) coastal current (not shown). At the mouth of the Lam river, we identified daily scenes where the outflow has the shape of a bulge with a recirculating current (not shown). There the Kelvin number is estimated around 0.35, so a recurrent bulge might possibly be formed. However, as in other delta systems, the typical length scale of the plume, and therefore the Kelvin number K, is not defined unambiguously, according to whether one considers the individual mouths or the system created by the merging plumes. Indeed, the plumes of the Red River distributaries (plus the Ma and Yen rivers) interact and merge in all clusters except cluster 3. Therefore, the overall fate of the Red River is likely to be better described by considering the system as a whole, as suggested by the 'menagerie' of plumes summarized by Horner-Devine et al. (2015; their Figure 5).

The Receiving Basin of Different Rivers in the GOT

As described in section "Observational Data Sets Used to Evaluate the Simulation, " the introduction of several tracers allows us to follow the pathway of different river systems: Red River, the rivers south of the RR delta (SR) and the rivers north of the RR delta (NR). This information is helpful for many practical purposes and environmental applications (e.g., contaminant dispersion, water quality). Mapping the 'receiving basin' of a river is a useful tool to estimate the zone of influence of dissolved or particulate terrestrial inputs, with the aim to design evaluation strategies of the coastal ocean ecological status and to build scenarios for the preservation or restoration of vulnerable areas. This is done, for instance, by Menesguen et al. (2018) in a numerical study to identify the source of eutrophication in the Bay of Biscay and the English Channel and to design a strategy of nutrient load reduction within European institutions directives framework. In this study, we define the "receiving basin" for different rivers from the passive tracer content with the same threshold as defined in section "River Plume in GOT: Identification and Variations of Area" (concentration > 7 unit/m 3 ), but this time we apply it to each river system. The calculation is made for the high runoff/summer monsoon period (July-August-September, hereafter JAS) and for the low runoff/winter monsoon season (December-January-February, hereafter DJF) from 2011 to 2016 (Figure 14). The choice to work with 'classical' seasons, and not on the clusters of section "Variability of the RR Plume, " is motivated by the concern to provide maps that can be easily used by the readers of this article. For each season, the "core of the receiving basin" is the area where the plume is present at least 90% of the time, the "average receiving basin" where the plume is present at least 50% of the time, "the envelope of the receiving basin" where the plume is present at least 10% of the time. Figure 14 also shows areas where the cores/averages/envelopes of the receiving basins for two different river systems superimpose.

In the low runoff period (DJF, Figures 14A-C), the receiving basins are very narrow and elongated along the coast. While the envelope of the RR receiving basin can extend to approximately 18.2 • N and mix with the plume of SR, its core appears mostly at the river mouths. In all cases, the receiving basins of NR is not detected with this choice of threshold for the tracer concentration (see below).

In the high runoff period (JAS, Figures 14D-F), the envelopes of the receiving basins of RR and SR are connected. South of 19 • N, the plumes of RR and SR are fully mixed and their envelopes are indistinguishable. North of 19 • N, the eastward extension of the RR envelope exceeds the one of the SR envelope and extends as eastward as 107.5 • E. The receiving basin of NR appears north of 21 • N but does not connect with the RR one. The average basins show different characteristics than the envelopes: while the SR influence can reach 18.3 • N, the RR basin does not extend further south than 19.3 • N.

In Figure 14, there is no connection between the RR and NR. However, it does not mean that they do not interact with each other. If we consider a lower threshold to identify the plume (e.g., 3 unit/m 3 ), the envelopes of NR and RR are connected both in low and high discharge periods (not shown). This suggests that the NR waters are strongly and quickly diluted. This also means that depending on the objectives of the study, and therefore on the choice of the tracer concentration threshold, the basins will be connected or not.

As written before, the information about receiving basins can help in case of water quality monitoring or contamination leakage. This analysis leads to 2 concluding remarks. Firstly, the freshwater along the Vietnamese coasts tends to flow southward and to exit the gulf. However, the freshwater from all the rivers tends to dilute to under 7% before it leaves the domain. This means that in case of contaminant leakage, the contaminants tend to mix and affect the gulf 's water quality locally mostly. Secondly, the RR average receiving basin extends further southward and in a wider coastal strip than the SR basin during the low runoff season, while the SR basin extends further southward and further eastward (south of 19.5 • N) than the RR basin during the high runoff period (Figures 14B,E). In both cases though, the RR runoff is higher than the SR runoff. These contrasted situations highlight the importance of wind and coastal circulation on shaping the receiving basin and consequently, on the fate of terrestrial inputs depending on their origin (SR or RR rivers systems).

CONCLUSION

We presented a comprehensive study of the plume formed by the waters from the Red River and three nearby rivers in the GOT using numerical simulation, in a realistic highresolution configuration over the period 2011-2016. Compared to various observational data, the model shows good results. We then compare several methods to identify the river plume in the study area. We found that identification through passive (dye-like) tracers is preferable since it allows to distinguish the runoff influence from the precipitation one and to distinguish the runoff from different rivers as well. The runoff shows large seasonal and interannual variability. It usually reaches the highest value in August and lowest in February, in phase with the monsoon system. At the surface, the plume area shows similar variations. However, the plume area reaches its peak in September, which is about 1 month later than the peak of runoff.

To identify the main spatial patterns of variability, we apply a clustering method to daily scenes of the model outputs where the plume is defined from the tracer concentration (>7 arbitrary unit/m 3 ). The cluster analysis identifies the plume regimes and their period of occurrence without having to pre-define these periods as one would do when computing seasonal averages for instance. Besides it allows identifying transition regimes (e.g., cluster 2).

The plume pattern of the first cluster is characteristic of the November-March period: the plume is narrow and contained within the inner shelf (bottom depth < 20 m), due to low runoff and downwelling wind (winter monsoon) with an intense downcoast current. The second pattern is mainly observed in the transition period of the monsoon (April, May, October) and is wider than the first pattern; both the relaxation of the winter monsoon and the weaker coastal current allow the plume to spread further offshore than in cluster 1. The third pattern occurs when the summer monsoon is the strongest and the runoff increases. The wind is upwelling favorable; the plumes are advected northward and detach locally from the coast. It is the only pattern where the plumes from the different rivers are disconnected. The final pattern coincides with the highest runoff period and the summer monsoon relaxation. The plumes are connected again in a pool of low salinity waters with the largest coverage, spreading both offshore and southward. Overall, we found that this hydrological system with multiple source points and a strong seasonal variability of along-shore wind is similar to the one for instance described by Kourafalou et al. (1996) in the Southern Atlantic Bight along the US coast.

The vertical variability of the plume is examined as well. In winter, the plume is usually mixed over the whole water depth. When the summer monsoon arrives at the end of March, the plume starts to detach from the bottom and spreads offshore as a surface buoyant layer. The plume can detach from the coast near 19 • N when a coastal upwelling develops. It may deepen offshore also. In particular, the plume is strongly deepened in the middle of August when trapped in a recurrent eddy near 19 • N.

The cluster analysis therefore evidenced at least three regimes. In winter regime, the plumes of all the Vietnamese rivers are connected in a well-mixed freshwater strip which is contained to the coast; the downcoast current is likely to act like a barrier preventing exchanges of riverine water and material cross-shore toward the center of the gulf. In early summer, the upwelling winds lead to a reverse circulation, where the plumes are advected northward and slightly offshore; the plumes are partly disconnected from each other. In the late summer regime, the wind relaxes and the large amount of freshwater due to the annual runoff peak spreads at the surface, reaching the maximal offshore and southward extension.

The impact of tides on the shape of the plume has been explored using a dedicated simulation without tides; we found that the general shape of the plume is not significantly influenced by tides, except in the Ha Long Bay area during the high runoff period (clusters 3 and 4). The main forcings of the farfield plume are the wind and basin-scale circulation. However, tides impact the vertical structure of the plume: they enhance the vertical mixing hence the riverine water concentration in the bottom layer. They also lead to a larger export offshore or the riverine water at the surface close to the mouth during the high discharge season. The influence of tides is obviously not limited to the far field; first, tides determine the buoyancy of the river water entering the ocean. Then, in the near and mid-field area, the bottom friction induced by tides is likely to impact the local dynamics and mixing, therefore influencing the plume properties in the far-field as well. An in-depth analysis of the tidal processes, of the tides-river flow interactions and of their impact on the velocity, shear and stratification in the estuaries and at the mouths was beyond the scope of this paper. It will be the topic of a future study aimed at better understanding the shape and dynamics of the outflow in the near and mid fields. The dynamics at the mouth is also highly important for sediment dynamics close to the mouth area.

Small-scale patterns of the coastal circulation, such as the recurrent eddy observed at 19 • N in August, are also involved in the dilution and fate of the riverine water in the ocean. Submesoscale features due to instabilities within the plume or at the fronts may also develop and influence the mixing as found in other systems (e.g., Horner-Devine et al., 2015;Ayouche et al., 2021 this issue). Such small-scale processes should be examined in dedicated studies.

At last, this study provides some reference information about the river plume variability, showing in particular the role of the wind. In such a context, ensemble simulations will be useful to assess the sensitivity of our results due to the uncertainty of the forcing conditions (wind and runoff).

Figure C1: total distance as a function of the number of clusters (line -linked to left axis) and amount of distance reduction (barlinked to right axis). The distance reduction

is calculated by subtracting the total distance of (i) cluster and (i+1) clusters.

To double check, we run KMA for two cases: 4 clusters and 5 clusters. The objective is to verify if using 5 clusters can lead to a new unique pattern, or if it just adds a pattern that is similar to one of the 4 existing patterns (obtained when we compute 4 clusters only). Let's call the patterns of the 4-cluster analysis A1, A2, A3, A4; and the patterns of 5-cluster analysis B1, B2, B3, B4, B5.

The results show that cluster A3 has similar pattern and temporal distribution as cluster B3. It is the same remark for clusters A1, A4, B1 and B4.

The difference is mainly visible in cluster A2. When we use 5 clusters, cluster B2 and cluster B5 are similar to cluster A2, but cluster B2 has a little bit less offshore extension and cluster B5 has a little more offshore extension than A2 (Figure C2). In short, using 5 clusters does not create a new different pattern compared to using 4 clusters. For all these above reasons, in this study, we decided to work with 4 clusters.

Cluster A2

Cluster B2 Cluster B5

Figure C2: cluster 2 in 4-clusters analysis (cluster A2), cluster 2 in 5-clusters analysis (cluster B2) and cluster 5 in 5-clusters analysis (cluster B5) Figure S5: frequency of occurrence of the plume within cluster 1 (A), cluster 2 (B), cluster 3 (C), cluster 4 (D) for GOT_REF (thick lines) and GOT_notide (thin lines) at 5%, 25% and 75%.

For instance, for cluster 1: at a given point, a frequency of 25 means that over the whole period when cluster 1 is present, the plume is present at this point 25% of the time.

IV.3 Further discussions

IV.3.1 The sensitivity of the clusters due to initial state and the reproducibility

In K-means clustering analysis, there are 3 problems that need to be addressed to have a reliable result. The first question relates to the number of clusters. This is addressed in the Appendix 3 in section IV.2. Since the K-means initial conditions consist of random numbers, the second problem is: "are the results sensitive to the choice of initial random numbers?", and the third question is: "if the initial condition is random, are the results reproducible or not?" In this section, I will address these 2 problems.

To produce results in section IV.2, in the initialization process, the first centroid of the cluster is taken "randomly" from the daily scenes of the plume area. Then, the clustering algorithm is run for 10 times, and the best result is kept. To check the robustness of the result, the process above is repeated 100 times. The results are shown in Figure IV.15. Most of the time, the clustering result is stable. A change of cluster is observed in less than 1% of the study period (15 days over 2192 days). Also, the change only happens 2-3 times over 100 times of run, with the maximum change being 5 times. Therefore, the result of the clustering analysis, in this study, can be considered robust and reproducible.

Figure IV.15

Highest frequency of occurrence (%). For example, 95% means that when the clustering analysis is run for 100 times (each time with 10 random conditions and the best result is kept), the plume pattern is attributed to the same cluster for 95 times, and to another cluster 5 times.

IV.3.2 The sensitivity of the tracer concentration near the boundary due to the sponge layer configuration

When looking at the surface tracer concentration on some days especially in the period of high runoff, it can be noticed that the tracer concentration shows a strong gradient near the boundary (Figure 4.16A). This strong gradient can be attributable to the sponge layer characteristics. Therefore, another twin simulation of GOT_REF, called GOT_S10, is created. The only difference between GOT_REF and GOT_S10 is that the sponge layer of GOT_REF contains 30 grid points, while GOT_S10 has 10 points.

Figure IV

.16 shows the surface concentration of GOT_REF and GOT_S10 on 15 August 2011.

In both cases, near the Southern boundary, the tracer concentration shows a sharp decrease to the South, but in GOT_S10, the area with low concentration is smaller. It means that the thickness of the sponge layer may affect the size of the low-concentration area near the boundary. However, it is not the only reason that creates this low tracer concentration area. In fact, it is due to the imperfection of the boundary condition and the dynamics of this region. In the model, the inflow tracer concentration is set to 0 at the open boundary. Therefore, if a tracer flows out of the model, it never gets replenished. In the summer, there is a recirculation at the south of the domain, as described in the Section 3 (Model assessment) of the paper above. A part of this recirculation is located outside of the model domain. Therefore, on one side, it brings water from the boundary (with tracer concentration = 0) into the model, but on the other side, it takes the tracer out of the domain. Therefore, the area of this recirculation is always filled with the water with tracer concentration = 0, which creates a very low concentration area near the southern boundary. This problem does not happen in the Hainan strait boundary because in the high runoff period the water direction is outflow.

Another test with the boundary condition is implemented. This time, the boundary condition of the tracer is equal to the concentration of the last point of the domain, and the other setting is the same as GOT_REF. However, this time, since all of the inflow water has the concentration of that last point, it leads to the accumulation of the tracer concentration inside the model and the results are more unrealistic.

III.4 Limitations due to the open boundary conditions

In GOT_NORIV (the simulation without any river), the same open-boundary forcing (same OGCM) is applied as in GOT_REF. Therefore, if there is an inflow from the boundary, and this inflow has a lower salinity than the salinity in the model, the low salinity water from the OGCM (created by river runoff in OGCM configuration) can flow into GOT_NORIV and create the uncertainty in the salinity of the GOT_NORIV.

To check if it happens or not, we looked at some sections at the coast near the southern boundary and found that the current is mostly southward, except between July and September, when an inflow can occur (Figure III.9). In the periods of inflow, we checked the salinity at some points along the section and found that the salinity from the OGCM is sometimes smaller than in GOT_REF. In these cases, some "unreal" freshwater can be indeed introduced into the domain so it can create uncertainty in the NORIV configuration. However, these differences are small (< 0.5 in all cases). This is the limitation for all the configurations without river if they use the same boundary condition as the reference configuration. To reduce this "unreal" freshwater, in the future, we can create another NORIV configuration with a larger domain, and use it as the boundary for the current NORIV configuration. Still, this can only reduce but cannot remove totally the unreal water flow into the GOT domain, since some unreal water can flow into the big NORIV. The salinity is taken at the circle for GOT_REF and at the square for OGCM. If the black line is positive and the colored line is negative, the unreal freshwater will be introduced to the model.

IV.4 Conclusion

In this chapter, we examine the variability of the Red River plume. The result has been published in the paper above. In summary, the river plume is first identified from a threshold of passive tracer concentration. This criteria is used to calculate the surface plume area. Then, we applied the K-means clustering method to daily scenes of the plume calculated above to identify the main plume spatial patterns. From sensitivity experiments, we have chosen to classify the plume patterns into 4 clusters, each linked with different forcing conditions. The vertical variability of the plume is examined as well, with different characteristics in different seasons. Finally, we tried to identify the "receiving basin" from different rivers.

The paper also contains maps of the first baroclinic Rossby radius of deformation for two seasons estimated from the model. To our knowledge this is the first regional estimate of this quantity at such a resolution in the GoT After the paper review, we did some further sensitivity experiments. Firstly, we tested the stability of the K-means clustering to the initial conditions.

The result shows that our results in the paper can be considered as robust. Then, we checked the boundary conditions and found that the largeness of the low concentration area near the boundary in the summer is linked with the parameterization of the width of the sponge layer. Finally, we examined the current and salinity near the boundary condition of the simulation without rivers and found that some low salinity water can be introduced into the domain in some day summers.

While the model's performance is almost good (compared to the data as in Chapter 3 and in the paper), these discussions above show that it still has some limitations at the boundary condition. This limitation can then lead to some uncertainties of the tracer concentration near the boundary. Beside the boundary conditions, there are other sources of errors. As shown in the clustering analysis, the plume pattern is partly controlled by the wind condition. Therefore, the errors in wind forcing may limit the conclusions of the paper regarding the plume regimes as identified by the clusters. One of the objectives of the next Chapter will be to study the impact of the wind forcing errors to the variability of the river plume. 2021), the wind plays an important role in the variability of the plume. Different wind regimes can create different plume variabilities, both in terms of surface coverage and thickness of the plume. A cluster analysis was performed in order to identify the main pattern of the Red River plume over the period 2011-2016. The authors found that in all 4 general patterns, the plume is strongly affected by the wind. Secondly, regarding the plume thickness, in winter, when the wind is downwelling favorable, the plume sticks to the coast and reaches the bottom. In summer, due to monsoon, the wind is upwelling favorable and even if the runoff is three times higher than in winter, the plume is thinner and advected offshore. show the area where the surface concentration is equal or higher than 7 unit/m 3 , which is the threshold to be identified as the river plume.

The work of [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF], as other studies on the river plume in other areas, rely on numerical modeling. However, the model results are known within some uncertainties. The uncertainties can come from the model itself (physical assumptions, parameterizations, imperfect bulk formula, …), from the model configuration (prescribed coastline and bathymetry, parameters, …) or from the forcing (atmospheric forcing, lateral boundary forcing, …). Understanding the uncertainties and their sources can help us better assess the model results in different scenarios, and quantifying the uncertainties is a necessary first step toward data assimilation. One way to estimate model uncertainties is to use one or several ensembles of wind uncertainty (amplitude, temporal and spatial structures). In contrast, [START_REF] Ayoub | Estimating uncertainties on a Gulf Stream mixed-layer heat budget from stochastic modeling[END_REF] estimated the error of the wind forcing by comparing it to a satellite product and then used the EOFs of the difference between the two products to create the perturbation. In this study, we will first attempt to estimate the amplitude of the uncertainty of the wind forcing by comparing the ECMWF product with a satellite product. In a second step, perturbations will be expressed and calculated as a pseudorandom linear combination of bivariate wind EOFs. Then, we will examine the response of the model to those uncertainties, at the surface and in the subsurface layer. We will revisit all the results that have been introduced in [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF], that is the surface plume area, the plume thickness and the stability of the plume pattern in terms of cluster classification.

We will also attempt to conduct an empirical assessment of the ensemble, that is to check if the estimated uncertainties from the ensemble are consistent with the model error deduced from model-observation misfits for the available observational data sets (as discussed for instance in Vervatis et al 2021a,b). To do so, we use rank histograms. This method allows us to assess the compatibility of the distribution of the ensemble with the distribution of the observational data.

The period of study is June to August 2015. It corresponds to the rainy season with high runoff and is a period when the wind goes through highly variable regimes. This chapter's outline is as follows. In section V.2, the method to create the wind perturbations is described. The results from the ensemble simulation are presented and examined in section V.3. The summary and conclusion will be delivered in section V.4.

V.2 Method

V.2.1 Numerical model and general configuration

The SYMPHONIE model (Marsaleix et al., 2006(Marsaleix et al., , 2008) ) is used in the same configuration (called GOT_REF) as described in Chapter II. The simulation runs from a save point on 04 April 2015 of the run which started originally from 01 January 2010. The perturbation starts from 01 May 2015 on. We then allow for a 1-month spin-up time. All the analyses are carried out from 01 June 2015 to 31 Aug 2015.

V.2.2 Ensemble generation

The ensemble is generated by perturbing the ECMWF wind forcing. The perturbations are meant to represent uncertainties on the ECMWF wind analyses. Those uncertainties being unknown, we need to make some assumptions to build the perturbations. Several methods are found in past studies, as reported by [START_REF] Ayoub | Estimating uncertainties on a Gulf Stream mixed-layer heat budget from stochastic modeling[END_REF] for instance. Here, we assume that the space and time structure of the wind uncertainties are the same as the space and time structure of the wind variability as in Le [START_REF] Hénaff | Assessment of observational networks with the Representer Matrix Spectra method-application to a 3D coastal model of the Bay of Biscay[END_REF]. This allows us to build the perturbations from the EOFs of the wind variability as described below.

Let's call W the 2D vector of wind velocity at 10m, composed of zonal U and meridional V components. We have:

𝑊(𝑥, 𝑦, 𝑡) = 𝑊 • (𝑥, 𝑦) + 𝑊 M (𝑥, 𝑦, 𝑡)

Eq. 5.1

With 𝑊 being the time-mean wind velocity and W' the variability of the wind.

We then form the Matrix of Samples (MoS) from the individual samples of W'. The wind EOFs are calculated by decomposing the MoS using SVD analysis. We end up with the EOF decomposition of the wind:

𝑊 M (𝑥, 𝑦, 𝑡) = t E *F+ 𝜙 * (𝑥, 𝑦) * 𝜃 * (𝑡)
Eq. 5.2

With:

𝜙 ' : bivariate (U,V) spatial component of EOF k (also called EOF k) -a nondimensional 2D vector field ("left" singular vectors) 𝜃 ' : temporal component (also called "amplitude") of EOF k -a dimensional 1D time series ("right" singular vectors)

The fraction of variance explained by EOF k, called Pk, is calculated as:

𝑃 8 = (𝑠 8 ) 9 ∑ : 8;< (𝑠 8 ) 9
Eq. 5.3 where 𝑠 ' is the k-th singular value.

The wind is perturbed by multiplying each of the summing operands in Eq. 5.2 by a factor (1+𝜎), where 𝜎 is a pseudorandom number drawn from a Gaussian distribution, with a zero mean and a pre-defined standard deviation 𝜀 (assumed independent of the rank of the EOF). This choice is made as we assume that the variance of the uncertainties represent about 𝜀 " *100% of the field variance.

The perturbed wind condition writes:

𝑊 % (𝑥, 𝑦, 𝑡) = 𝑊(𝑥, 𝑦, 𝑡) + ] E 'D.
𝜙 ' (𝑥, 𝑦) * 𝜃 ' (𝑡) * 𝜎 𝒌 (𝒕) Eq. 5.4

Values of 𝜎 are pseudo-randomly drawn for each EOF rank and each member, and are updated every 5 days, which is equal to the value used in the Bay of Biscay in the previous studies (Quattrochi et al. 2014[START_REF] Ghantous | Ensemble downscaling of a regional ocean model[END_REF]; in between, they are linearly interpolated. In this study, we also update the pseudorandom numbers every 5 days, as it reflects well the synoptic time-scale over the region (e.g. tropical cyclones, cold surges).

The perturbations are computed at the time step of the ECMWF fields, that is every 3 hours. They are then linearly interpolated at every model time step and spatially interpolated onto the model grid.

Another parameter of the ensemble simulation is the number of members. We ran our simulations with a total of 50 members, divided into 25 pairs. In each pair, the absolute value of the perturbations are the same, but the signs are opposite. Therefore, at any time, the ensemble average of perturbations is 0 (but the ocean response can have nonzero ensemble mean due to nonlinearities in the model).

To sum up, to create the wind perturbations, we will need the wind EOFs and pseudorandom numbers. In the next 2 sections, we will present the choices we will apply for both of them .

V.2.3 EOF wind analysis

Figure V.2 shows the spatial and temporal components of the first 5 bivariate (U,V) EOFs and the fraction of variance explained by the first 12 EOFs, calculated from 3-hour ECMWF wind from 01 May to 31 August 2015 over the GOT domain. The EOF decomposition is computed using SVD (see Eq. 5.2) as implemented in the NumPy library2 in Python. From 19 -26 June, there is a tropical cyclone in the gulf. It is an extreme event and not in the scope of this study, therefore we decided not to use winds in that period to calculate the EOFs. It means that in that period, we chose to set the perturbation to 0 (therefore, in the ensemble simulations, the model will run with unperturbed wind forcing). (C,D),(E,F),(G,H),(I,J): same but for EOF 2, 3, 4, and 5, respectively. (K): Fraction of variability explained by the first 12 EOFs. The upper number is the ratio of the EOF variance to the total variance, while the lower number in the bracket is the cumulative variance.

EOF 1 is the main contribution to the wind, explaining 49.3% of the total variance. Its spatial component shows a wind mostly in the N-S direction, following the shape of the gulf. The velocity is lower near the coast and reaches its highest values near the center -north of the gulf. The temporal components in May and June show negative values, which means that the wind direction is mostly southerly. From July, the wind direction reverses with positive values being more observed. EOF 2 accounts for 21% of the total variance, representing mostly a zonal wind pattern with the highest velocity located at the north of Hainan Island. EOF 3 explains 5.1% of the total variance, with weakest values at the center of the gulf (near 20°N), and highest at the coast.

The wind perturbations are built from Eq. (5.4) with a limited number of EOFs. Here, the first 5 EOFs are kept to reach at least 80% off the total variance. After EOF 5, the percentage of accumulated variance does not change much (from EOF 5 to EOF 12 it only increases from 80.3% to 88.6%). EOF 4 shows a very strong wind near Hainan Strait. We suspect that this pattern may not realistic, therefore, in this study, 4 EOFs are used (1,2,3,5) with total variance ~ 77.5%.

V.2.4 Random error generation

Another challenge of the perturbation generation is the choice of the amplitude of random errors. Quattrocchi et al. ( 2014), Vervatis et al. (2016) and most recently, [START_REF] Ghantous | Ensemble downscaling of a regional ocean model[END_REF], in their stochastic studies in the Bay of Biscay, create the ensemble with a Gaussian pseudorandom number generator with standard deviation 𝜀 = 0.3, which assumes that the standard deviation of uncertainty on the wind velocity is 30% of the standard deviation of the wind. However, this value is not necessarily general as it is based on the assumptions of the authors for a specific situation. [START_REF] Pasmans | Ensemble of 4DVARs (En4DVar) data assimilation in a coastal ocean circulation model, Part I: Methodology and ensemble statistics[END_REF] estimate 𝜀 from satellite and in situ observations of wind and from a Bayesian Hierarchical Model; in their case, 𝜀 differs according to the rank of the EOFs. They found that in their area of study, the value 30% is usually an underestimation, especially for high modes. In this study, we propose to estimate 𝜀, using ASCAT satellite data as well but in a much simpler approach and without introducing modal dependency.

The ASCAT product distributed by EUMETSAT (EUMETSAT, OSI SAF, 2021) contains level 2 wind data, retrieved from the Advanced Scatterometer (ASCAT) on MetOp-B satellite, at 12.5 km sampling resolution (with an effective resolution of 25 km). It is available daily and can be downloaded freely from the JPL website: https://podaac.jpl.nasa.gov/dataset/ASCATB-L2-Coastal . The data are downloaded in the period from 01 May to 31 August 2015. The ASCAT User Manual indicates that 'the accuracy is better than 2 ms -1 in wind component standard deviation with a bias of less than 0.5 ms -1 in wind speed'.

Before doing the calculation, the data are processed through several steps. Firstly, since the data have different data tracks on different days, they are interpolated onto the ECMWF grid. Then, the data that have less than 50% spatial coverage are filtered out; over the 123 days of this period, only 72 days (~58.5%) are used for further calculation. Finally, all the missing values of the ASCAT are mirrored to ECMWF data values to make the two datasets homogenous.

Then, the spatially-varying daily mean bias between the ECMWF and ASCAT is calculated, and the temporal standard deviation (hereafter std) of that bias is shown in Figure V.3A (for U component) and Figure V.3D (for V component). For both U and V, the std of the bias is higher near the Red River coast, reaching around 3ms -1 for U and 3.5 ms -1 for V (note that these values are larger than the accuracy of the ASCAT product). The temporal standard deviation of the ECMWF wind is shown in Figure V.3B (for U) and Figure V.3E (for V). As for the bias, the maximum variance of ECMWF V is higher than U. However, the maximum variability in U is shown at the south of the gulf (~3.7 ms -1 ), while for V it is at the north of the gulf (~ 5 ms -1 ). The ratio between the standard deviation of bias and standard deviation of ECMWF is shown in Figure V.3C and V.3F. For both U and V, this ratio is higher at the Red River coast and can reach a maximum local value higher than 1. On average, this ratio (a proxy for 𝜀) is around 0.8 for U and 0.7 for V.

Another way (not shown) of calculating the standard deviation has been considered also. Instead of calculating the ratio between the two maps of std(bias) and std(ECMWF), then taking the spatial mean, we just calculate a single value for the std(bias) and std(ECMWF) over time and space and then take their ratio. This time, the ratio is 0.8, 0.6, and 0.6, with respect to U, V and using simultaneously U and V.

In conclusion, these comparisons show that it seems reasonable to choose the value of 𝜀 within the range [0.6-0.9]. However, (1) if the perturbation is too large, the perturbed wind can go too far from the reference dataset and may totally change the reference pattern, and (2) we did not take into account the ASCAT measurement error in calculating the ratios. As a compromise, we choose 𝜀= 0.6. To further assess the appropriateness of this choice, we examine what we could have obtained with another method. In the absence of information from observations, we can use wind fields from another year or from a typical year. In that case, the perturbations do not mimic the wind uncertainties; instead, the sensitivity of the ocean solution is assessed with respect to other 'wind solutions' (in the space of the ECMWF model's solutions).

Here we compute the daily mean wind fields from 10 years of ECMWF wind field (hereafter ECMWF_MEAN). The bias is computed between the year 2015 and ECMWF_MEAN and then we do the similar calculation as we do with ASCAT, as if ECMWF_MEAN were observational data. In this case, the interannual variability of the wind is considered as the uncertainty of the wind. Such an approach is often adopted in climate sciences: to evaluate the skill of a prediction scheme, the prediction is often compared to the climatological value in order to evaluate if the forecast gives some added values compared to the simple climatological mean. The standard deviation of the bias for V is higher than U, ranging from 3 ms -1 at the south of the gulf to 5 ms -1 at the north. The ratio between std(bias) and std (ECMWF) varies slightly, ranging from 0.9 to 1.0 in almost all areas. These calculations confirm that the choice of ε should be higher than 0.3, the value chosen in other studies. For the reason given earlier, we take ε =0.6. In the following we however present some results from an ensemble generated with ε =0.3 as a reference.

V.2.5 List of simulations

The simulation without wind perturbation is called GOT_REF. There are two ensembles. The first ensemble, called ENS_REF, which uses 𝜀 = 0.6. Another ensemble, called ENS_0.3, has the same configuration as ENS_REF, but with 𝜀 = 0.3. The wind perturbations are computed online during the SYMPHONIE simulation. -2 , except the first half of June and on 01 August. For the V component, the unperturbed wind has two significant peaks: the first peak is the storm mentioned before, from 23 -25 June. In this period, the maximum wind stress can reach 0.17 Nm -2 . The second peak is observed from 20 July to 31 July. In this period, the maximum wind stress can reach around 0.11 Nm -2 , which is lower than during the storm, but this event lasts longer than the storm. For both components, the spatial mean of ensemble spread is less than 0.03 Nm -2 most of the time, except for the period of 05-11 July. The peak of the spread of wind stress in this period is consistent with the peak of EOF 1, suggesting that the spread is sensitive to the particular pattern of EOF 1. However, in another peak of EOF 1, such as from 21 -25 August, the correspondence is weaker, suggesting more complex effects. The spatial pattern of ensemble spread on July 10 shows a north-south gradient as well (Figure V.6C) but stronger, with the smallest spread value at the south (around 1 cm) and largest at the north (30 cm). This coincides in time with the maximum of EOF 1 (wind perturbations mostly along the N-S direction, translating into possible water piling up/entrainment). Interestingly, it has been known from previous studies (eg Nguyen et al., 2014) that the diurnal tides are resonant which leads to their amplification in the northern GOT. We could assume that a large wind stress perturbation with a time scale of about 1 day (such as those on July 10) would generate a perturbation of the sea level that would be resonant and amplified with a similar structure as the diurnal tides (see for instance The ensemble spread for 09 July is shown in Figure V.8 A, B. On this day, the spread of U is highest near the RR coast (more than 20 cm s -1 ), and the spread of V is highest near the west of Hainan Island (more than 30 cm s -1 ). The high spread along the coast in V component suggests that the wind perturbation has an impact on the coastal current created by the river runoff.

Near the RR mouths, the maximum value of spread is more than 20 cm s -1 for both components ( Near Ma River, the river configuration is the same, however since the spread of the surface current is very small (Figure 8C, D) therefore the spread of the SST near this river mouth is also small.

Figure V.12D

shows the spread near the RR mouths, where the spread is the highest. At the river mouth, the spread is very low due to the homogeneity of the riverine water properties which are set by the runoff and the tides. Downstream, the spread increases and can reach 2 near the Yen river. The high spread near the rivers can be due to 2 reasons. Firstly, wind perturbation can change the wind-driven current at the surface, and then change the horizontal advection. Secondly, it can have an impact on the vertical mixing. In that case, the uncertainties on the plume should also be significant at the subsurface. In the next section, the subsurface spread will be examined.

V.3.1.2 Ensemble spread at the sub-surface

The vertical section at 19°N of the ensemble spread for temperature on 11 August is shown in and [START_REF] Ayoub | Estimating uncertainties on a Gulf Stream mixed-layer heat budget from stochastic modeling[END_REF]. The change of wind can have both impacts: at the surface and at the sub surface. However, since the horizontal gradient of the temperature is small, even if there are some horizontal displacements, they do not have any signature on the change of the temperature.

For salinity, the maximum spread occurs near the surface (Figure V.13D), where the salinity is low due to the runoff. Similar to the temperature, the maximum spread should occur where the gradient is highest. For salinity, the horizontal gradient is high where the river plume is present, therefore the spread is highest at the surface. The vertical gradient of salinity is small except for the area strongly influenced by the runoff. Here, the condition is reversed to the temperature. Even if the wind perturbation changes the vertical mixing, the displacement of the water mass with similar salinity does not create any signature on the salinity ensemble spread.

The vertical distribution of the spread for the salinity and temperature shows that the perturbations of the wind can have an impact both at the surface and subsurface. Depending on the variables, the spread may be higher at the surface or at the subsurface. On average, the spread of the plume area is around 1452 km 2 , which is equal to 18.8% of the plume area in GOT_REF. The peak of plume spread is from 20 -20 August at more than 3600 km 2 , but the percentage of the spread reached the maximum value on 02 June (more than 50%). The relationship between the spread and the GOT_REF plume area is not clear, except in the end of August where both the spread and the GOT_REF plume area reach their peaks. Interestingly, the members have a different behavior from August 13-15: in some members the plume area continues increasing, while in others, it decreases more or less (as in GOT_REF) and, in some cases, stays below their August 15 value until the end of the run.

Similarly to the spread of the salinity and temperature, in general, the high spread of the plume can be due to the change of the vertical mixing and surface current.

V.3.2.2 Clustering analysis

As in Chapter IV, to assess the uncertainty of the plume pattern, we classify the daily plume pattern of each member of the ensemble into 4 clusters, based on the reference pattern classification in [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF]. In other words, we use the results from 6-year cluster analysis as the training data for this study.

To do that, the daily plume pattern of each member is identified by the threshold as described in Chapter IV (>=7 unit/m 3 ). Then, the distances between this daily pattern and the centroids of the 4 reference clusters described in [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF] are computed. Finally, the daily pattern for each member of the ensemble is attributed to the cluster which has the minimum distance.

We calculated a simple Maximum Likelihood estimate by choosing, for each date, the cluster that occurred the most frequently across the ensemble (denoted as MAX estimate). 2021) are yet unchallenged and robust with respect to wind perturbations.

In GOT_REF, in this summer period, cluster 3 occurs the most (54/92 days) and also is the most stable cluster, in the following sense: across the ensemble, we do not observe many changes from cluster 3 to other clusters; if it changes, it mostly changes to cluster 1 (5 July) or cluster 2 (5 August).

Cluster 1 happens in just 2 days (06 and 07 July), and fluctuates across the ensemble between clusters 1 and 2. In [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF], cluster 1 is described as the pattern that occurs in the lowest runoff period, so it is not expected to occur much in the rainy season in July. Its fluctuation between cluster 1 and 2 when submitted to the wind perturbations suggests that the plume pattern on 06 and 07 July in GOT_REF is on the edge of cluster 1 and cluster 2 and just closer to cluster 1. In that case, the attribution to cluster 1 may not be as robust as the one to cluster 3.

The plume pattern appears to be the least robust to the wind perturbations in 3 periods: 10-22 July, 5-12 August, and 21-28 August. It is the period where we can see at least 3 different clusters occurring in different members on the same day. The high variability of the plume patterns in those periods is consistent with the large spread of the sea surface salinity, suggesting that the runoff could be another major source of uncertainty along with the wind. We will come back to this in the conclusions. Though Figure V.15 shows that most of the time, amongst 50 members, the most frequent cluster classification is the reference cluster classification, it also shows that sometimes the wind perturbation can significantly change the regime of the plume. For example, in member 33, the plume is classified as cluster 3 most of the time, though the reference is cluster 1 (e.g. 6-7 July), cluster 2 (e.g. 8 -15 July), cluster 4 (e.g. 16 -21 July).

To confirm that the wind perturbation can indeed change the plume regime, and clarify what cluster 3 is able to capture in member 33, here we analyze two short periods.

From 5 -10 July, member 33 is labeled as cluster 3. The daily mean surface concentration shows that the plume is advected northward, and the plume from Ma river is disjoint from the plume from other rivers, which is indeed the typical characteristic of reference cluster 3 (Figure V.16). North of 20N, the wind is stronger and favors a northward advection, which explains the offshore transport of the plume. South of 20N, the wind is weaker and the direction changes to northerly, however, the current is northeastward so the plume from Ma river is advected offshore. In member 34, which belongs to cluster 1, the plume is thin and attached to the coast, which is also the typical shape of cluster 1. The strong northeasterly and westerly winds in the north and south of the GOT respectively trigger the southward coastal current, which explains why the tracer is trapped near the coast.

The similar situation happens again from 15 -20 July (Figure V.17). In member 33, the surface plume is advected northward, consistent with wind forcing and Ekman transport, while in member 34, the coastal wind is weak and the current is southward, which explains the thin plume attached to the coast. The shape of the plume in member 34 is similar to the typical shape of cluster 2. To better estimate the spread of the plume thickness, we look at the points at which the plume is present in at least 30 members (dotted red in Figure V.18B). In June, the spread is less than 1m. The peak of the spread happens from 09-13 July (> 3m), which is the same time of the largest spread of wind stress, suggesting that the change of the wind condition also has a considerable impact on mixing. Other high values of the spread occuralso in the period of 5-12 August and 21-28 August at the same times as the large variations of the plume pattern (Figure In the period of 20-25 July, the plume in GOT_REF reaches the highest thickness, increasing from 3m near the coast to 10m offshore. The deepening of the plume is associated with an eddy (Figure V.19A), such as the one described in [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF]. The map of relative surface vorticity on 24 July 2015 shows the signature of an anticyclonic eddy, and the plume is deepening inside. In this period, the spread is around 1-2m. Our first idea is that the processes causing or involved in the eddy have a critical impact on the mixing of the plume and the wind perturbations do not change much the mixing inside the plume. However, the plume thickness computed from each member of the ENS_REF shows that the wind perturbation can alter the current, change the shape of the eddy and finally modify the plume thickness. Another interesting thing is that the eddy in GOT_REF is located in the area where the plume is present in less than 20 members of ENS_REF, which means it is unstable.

V.3.3 Event analysis

In this section, we will examine the impact of the wind perturbation on the plume variability in some special events. Along the coast from 20 -21°N, the plume seems to be stable across the ensemble, as demonstrated by the small difference between the 5-member and 45-member plume occurrence. Along the latitude of 20N, the plume spread is large. For half of the members, the plume can extend to 106.3°E, while in 5 members it can reach around 107.4°E, which is more than 100km from the coast. A large spread is also observed near 19°N.

Figure V.21 A, B, C

shows the surface tracer concentration of GOT_REF, of the member with the largest plume area (member 41) and with the smallest plume area (member 42) on 27 July. In all 3 simulations, the passive tracer tends to extend offshore, however, the surface tracer concentration of member 41 is the highest and member 42 is the lowest.

The daily wind condition corresponding to each simulation is shown in Firstly, it changes the mixing. A stronger wind increases the vertical mixing, therefore the tracer at the surface is decreased and the tracer concentration at the sub-surface is increased.. Indeed, the histogram of vertical diffusivity shows that the mixing is stronger under the stronger wind Secondly, the wind perturbation also changes the surface circulation. Near the RR coast, the surface current is eastward and northeastward in all cases (Figure V.21 M, N, O). However, the current is the strongest in member 42 in which the wind is the strongest. Under stronger current, the plume passive tracer is advected further and surface concentration is more diluted. The concentration of tracer can eventually fall below the threshold of the plume.

V.3.3.2 20 August

From 15 -20 August, the river plume shows a large spread (Figure V.14A), although the spread of the wind stress is small (Figure V.5A). However, it is the period where the runoff is high. It means that the runoff plays an important role on the spread of the plume. On July 27, when the runoff is low, the increase of the current due to the wind perturbation can dilute the plume, leading to a smaller plume area. However, on August 20, when the runoff is high, the increase of the current can help the plume to be advected more offshore and increase the plume area.

Using only the wind condition on August 20, it is hard to assess the effect of wind on mixing.

In the first 15m, member 39, which has the strongest wind, seems to have the strongest mixing, as shown by the vertical profile of concentration at 20°N (Figure V.23L). However, from 30m and below, the tracer concentration of member 40 is higher (Figure V.23K). It means that the dynamic on a day does not only depend on the current conditions, but also on the past conditions, which have been changed due to previous wind perturbations. Member 39 and 40 have been generated by paired, opposite-sign perturbations, but have totally different surface currents and tracer concentrations. The changes do not happen in just one day. Therefore, in the future, the impact of the wind perturbation could be assessed with enhanced benefit over a longer period. To assess the ensemble with respect to observations, rank histograms (hereafter RH) are used. This method allows us to assess the compatibility between the distribution of the ensemble and the distribution of the observational data. It is also called the Talagrand diagram.

The principle of this method is as follows. In the ensemble simulation, the probability density function (pdf) of a variable X is estimated from N values Xi (i = 1,...,N) with N being the number of members of the ensemble. Ranking these values in ascending order defines N+1 intervals (ranks). The observational data Yk (k = 1,...,P) are then placed on the suitable rank based on their values. We use the following convention: if Yk is smaller than all ensemble members, it is assigned as rank 1, reversely, if it is larger than all ensemble members, it is placed at rank N+1. The histogram of the position of Yk with respect to the Xi 's defines a measure of the statistical consistency of the ensemble with respect to the observations: If Yk is an independent realization of the same pdf which has produced the Xi 's, it will be statistically indistinguishable from the Xi 's, and will therefore fall with equal frequency 1/(N + 1) in each of the intervals defined by the Xi 's [START_REF] Talagrand | Evaluation of probabilistic prediction systems[END_REF].

The flat rank histogram indicates that the distribution of the observation is well represented by the ensemble (Figure V.24). Other than that, a rank histogram with a U-shape indicates that the spread of the ensemble is not large enough (the ensemble is underdispersive), while a dome shape means that the ensemble spread is too large (overdispersive ensemble). The figure also shows biased ensembles. Note that types can be mixed, e.g. underdispersive and high-biased.

Another remark is that the presence of observational errors will modify the observational pdf and can make conclusions misleading. Rank histograms are often used in ensemble/stochastic modeling and/or data assimilation. It can be used to empirically assess the ensemble, e.g. in [START_REF] Vervatis | Assessment of a regional physical-biogeochemical stochastic ocean model. Part 2: Empirical consistency[END_REF] and [START_REF] Yan | Ensemble assimilation of ARGO temperature profile, sea surface temperature, and altimetric satellite data into an eddy permitting primitive equation model of the North Atlantic Ocean: ENSEMBLE ASSIMILATION TO AN OCEAN MODEL[END_REF] with respect to the physical (SST, SSH) and biological variables (phytoplankton functional types, chlorophyll). [START_REF] Santana-Falcón | Assimilation of chlorophyll data into a stochastic ensemble simulation for the North Atlantic Ocean[END_REF] use this method to assess the performance of the free-run and assimilated ensembles with respect to chlorophyll observations from SeaWiFS on selected regions. [START_REF] Van Velzen | OpenDA-NEMO framework for ocean data assimilation[END_REF] use rank histograms to assess the performance of the Ensemble Kalman Filter (EnKF) data assimilation on SSH in different time scales and different numbers of members of the ensembles.

Here, we compute the rank histograms for SST and for surface currents, focusing on the plume area. In accordance with the usual practice for the calculation of rank histograms, the observations are not perturbed.

To compute the rank histograms of SST from ENS_REF and OSTIA data, firstly, for each day in June-August, the ensemble SST is spatially interpolated onto the OSTIA grid. Then, at each point of the area between 105.5 -107.5E, 18.5N -21.5N (Figure V.25A), OSTIA data are attributed to the suitable rank according to the procedure indicated above. The results are shown in Figure V.26. Often in the period, the histograms appears U shaped, with more values at rank 1 and rank 51, suggesting that the ensemble is underdispersive when this happens. However, in generating the ensemble, we only took into account one error process (surface wind), while it is not the only error process at play, on one hand, and on the other hand the OSTIA observations may contain errors.

The histograms typology is found to vary with time. In June, the ensemble is quite consistent with the OSTIA data, as shown with the relative flatness of the histograms, except on 3 days at the end of the month when it is low-biased (many data project onto class 51). The ensemble is the most low-biased from around July 25 to August 08. This inconsistency may be not related to the wind uncertainty but to the limitations of the river runoff configuration in the model. In general, the ensemble is found to be underdispersive, although results vary with time. For the U component, it is U-shaped in June, but it is high-biased in mid July and low-biased in middle and late August. For the V component, the ensemble appears low-biased in the first half of June, then the histograms have a U shape most of the time.

For both SST and surface currents, the histograms flatness seems to be more pronounced in July 5 -15, which is the same time as the large spread of wind stress. This could mean that increasing the perturbation may improve the performance of the ensemble at describing the observed variance in statistical terms. The y-axis shows the rank. There is no data on 20 and 21 June.

V.3.4.2 Uncertainties budget

Another way to assess the relevance of the ensemble in simulating the model errors is to evaluate the skill of the spread at reproducing model errors. This method is based on the fact that we have two estimates for the model errors: the ensemble spread and, on the other hand, the misfits between the model and the observations. Of course the comparison is only valid within the observational errors. This test is useful in particular to assess the uncertainty budget in ensemble-based data assimilation methods.

In concrete terms, we first compute the spatial mean of the ensemble spread for each (U, V) component over the HF radar domain. We also compute the difference at each radar grid point between the ensemble-mean current and the observed one (for each component); we then compute the RMS of these differences over the radar domain (hereafter 'the RMSE').

Figure V.27 shows the temporal evolution of the spatial mean of ensemble spread and the ensemble-mean RMSE with respect to the HF radar current data. For the U component, most of the time, the RMSE is around two times larger than the ensemble spread. On average, the temporal mean of RMSE and ensemble spread are 8.05 cm/s and 4.08 cm/s, respectively. The ensemble spread varies from 2 -6 cm/s and increases gradually from June to August. The RMSE ranges from 4 -15 cm/s (except the period of the tropical storm), with a peak in late July and early August. However, the difference between the two curves stays well within the observational error estimates: considering the end July/beginning of August period, taking 25cm 2 /s 2 for the ensemble variance in U and 100cm 2 /s 2 for the MSE of U, we end up with a variance budget of observational errors of 100-25=75cm 2 /s 2 , i.e. a st.dev. of 8.66cm/s, which is compatible with the HFR observational error estimates. Indeed, Tran et al. (2021) estimate that the RMSE of HFR data with respect to the current measured at a mooring during 14 days in October 2015 is 12 cm/s for zonal current and 10 cm/s for meridional current (as shown in Chapter IV).

For the V component, the mismatch between the RMSE and ensemble spread is smaller. On average, the temporal mean of RMSE and ensemble spread are 7.16 cm/s and 4.76 cm/s, respectively. Similarly to the surface current rank histograms (Fig. V.26), the ensemble spread is increased and well matched with the RMSE from 5 -15 July, which is the period of large wind stress spread, consistently with the result that the current is sensitive to the large wind perturbations. As for the U component, the difference between both curves is compatible with the estimated HFR observational errors, although the curves are closer to each other this time, and therefore the resulting observational error variance budget is smaller. This is in fact consistent with the smaller RMSE of the V component as estimated by Tran et al. (2021).

V.4 Summaries and conclusions

In this Chapter, ensemble simulations are used to assess the impact of an important source of uncertainties, namely the wind uncertainty, on standard model variables in general and on the variations of the Red River plume in particular. In particular, we aim at assessing the robustness of the Chapter IV results and conclusions with respect to wind forcing uncertainties. In addition, we tried to address the following more general questions with the ensemble approach:

1. what is the space and time variability of the ensemble spread of the surface ocean variables ?

2. can we identify the physical processes at work to explain the sensitivity ?

One preliminary methodological question was "how to generate wind perturbations that are indeed representative of ECMWF wind errors ?". By comparing the daily ECMWF data to the ASCAT data and ECMWF climatological data, the uncertainty of the ECMWF wind in this area is estimated as 60% of the wind variability. The wind perturbations are then created as linear combinations of bivariate wind EOFs weighted using a Gaussian pseudorandom number generator with standard deviation equal to the uncertainty value estimated above. The ensemble is run with 50 members, perturbed by pairs of opposite perturbations so that the ensemble average of perturbations is zero.

The ECMWF to ASCAT comparison also indirectly answers another methodological question, "is the generated ocean ensemble significantly representative of model errors ?", insofar as we calibrated in this manner the stochastic forcing. But the question also needs of course to be addressed for the ocean variables (see below).

Another preliminary methodological question was "how (by means of which metrics) do we define the sensitivity to wind uncertainties with an ensemble? In terms of metric definition?".

In the adopted approach, the analysis of the ensembles in terms of sensitivity consisted mainly in calculating the ensemble standard deviation (the "spread"), as the ocean response of ocean variables to the stochastic forcing, and the statistical distributions of ensemble samples.

The assessments are mostly based on the ensemble spread, which is computed as the standard deviation along the ensemble with respect to the ensemble mean. First, we assess the spread of the model variables. The results show that the SSH spread is less than 2 cm most of the time.

The spatial mean spread of the SSH over the whole period decreases southerly. The spread increases locally near the coast where the river plume is present. The peak of the spread happens at the same time as the peaks of the wind stress, which is also the peak of the EOF1, indicating that the result may be sensitive to the EOF analysis.

The surface current shows the peaks of the ensemble spread mostly at the same time as the peaks of SSH spread, although the links between both would have to be further clarified. The high spread along the coast in V component suggests that the wind perturbation has an impact on the coastal current created by the river runoff.

The temporal evolution of the spatial mean spread of SST shows the same time pattern as the SSH. However, unlike the SSH, the maximum spread of SST is highest near the south of Hainan island. We expect that the spread of the wind leads to the change of the mixing (which is already high due to tidal effect in that area) and then leading to the displacement of the water mass that has a high temperature gradient. Over the whole period, the mean spread is highest over the Red River coast, which can be related to the difference between temperature of the riverine water or the upwelling present in this area. The SSS, as expected, has high spread near the Red River.

The vertical section at 19°N is investigated to assess the vertical structure of the spread of temperature and salinity. The temperature spread is highest around 20-30m deep, where the vertical gradient is highest (base of the mixed layer; as in [START_REF] Andreu-Burillo | Assessing oceanmodel sensitivity to wind forcing uncertainties[END_REF][START_REF] Ayoub | Estimating uncertainties on a Gulf Stream mixed-layer heat budget from stochastic modeling[END_REF], a clear confirmation of the fact that the uncertainty of the wind affects vertical mixing. On salinity, the spread is highest near the surface in the vicinity of the sanity gradient due to the runoff.

Next, the spread of the Red River plume is assessed. At the surface, the spread of the plume area is highest in August, but it does not follow the same trend as the SSS. The maximum spread of the plume area happens at the same time as the maximum of the plume area, suggesting that our results might be sensitive to the uncertainty of the river runoff as well.

The clustering analysis was carried out per-member using as a reference the clusters of GOT_REF, as computed in Nguyen Duy et al. (2021). Results evidence some changes of clusters across the ensemble, but almost all the time, the pattern that is most likely to occur is still the reference cluster. As a consequence, we estimate that within our experimental conditions the findings of Nguyen-Duy et al. ( 2021) are yet unchallenged and robust with respect to wind perturbations.

Finally, the spread of the plume thickness shows that the plume is not only changed vertically but also changed horizontally. In June, the spread is less than 1m. The peak of the spread happens from 09-13 July (> 3m), which is the same time as for the largest spread of wind stress, suggesting that the change of the wind condition also has a considerable impact on mixing. The plume can deepen inside the eddy, as discussed in [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF]. The wind perturbation can alter the current, change the shape of the eddy and finally modify the plume thickness.

The spread of the plume for some specific days is also examined. The stronger wind increases the vertical mixing. In particular, stronger upwelling favorable wind can increase the current and help the plume to be advected further. In two members which have in-pair perturbations, the current conditions can be totally different in both direction and value, which is due to the past conditions, especially in GOT since the wind can change very fast. Therefore the impact of the wind perturbation should be assessed over a period (not instantaneously).

We would like to finish with several remarks and limitations of our approach.

First, we found out that the wind difference patterns obtained when matching ECMWF data with ASCAT data were not quite the same as the patterns generated by our base of EOFs, even if the error scales seemed similar. We did not take into account the spatial structure of ASCAT errors when doing that comparison, because we do not know those errors. As a consequence, we limited ourselves in this study to using that comparison to get a very rough estimate of ECMWF error variances, at the risk of overestimating them and of structures of smaller scales than those actually taken into account, and did not try to infer anything from the difference patterns themselves.

Second, in this study, we limited the source of stochastic perturbations to the surface wind only.

Although several studies, cited in the introduction of this Chapter, agree on the wind being an essential, if not dominant, source of ocean uncertainties, we could have considered other sources. For instance, we became aware during this study that the runoff could be another major source of uncertainty impacting the RR plume, along with the wind. As another example, we could not fully clarify the apparent partial correlation between the spread of surface currents and the SSH spread -in order to better explore that link, atmospheric pressure would need to be perturbed in a consistent way with the wind, in particular since the GoT is a shelf sea and semi-enclosed sea, so it should be sensitive to water displacements as, e.g., a low-pressure center passes through. One way to do this comparison in a meaningful way would be to use multivariate ECMWF ensembles, perturb both the wind and the pressure, and perhaps add a wave model to better model the wind effects on currents. We did not have the time to explore those two directions, but it will need to be done.

4. We ran an ensemble of simulations assessing the model response due to perturbations of the wind forcing. We analyzed the variability of the resulting uncertainties over space and time of the ocean variables. The K-means clustering approach was used on ensemble members, providing an assessment of the robustness of the results of 3. above.

Using the tool described above, combined with the observational data, we obtained the main results below:

Model assessment

The model results from GOT_REF have been assessed by comparing with available observational data: satellite data (altimetry for tide, OSTIA for sea surface temperature, SMOS for sea surface salinity), in situ data (VITEL, CFO, both datasets with temperature and salinity profiles), and HF radar (for surface current). We also compared our results with the OGCM (the global circulation model) to emphasize the added value of our high resolution regional model compared to a global model. Overall, despite some biases, the model reproduced well the seasonal conditions even in the very shallow area. The comparison with the CTD data near the coast shows an overestimation of the surface salinity, therefore we suspect that the runoff may be underestimated.

We also discussed the typical horizontal length scale (first baroclinic Rossby radius) computed from the simulated density field and described the general circulation.

RR plume identification

To identify the RR plume, different methods were taken into consideration. Firstly, we applied a method based on the Stratification Index (hereafter SI), assuming that the river water forms a buoyant layer over the ocean water. The second method consisted in identifying the plume from the sea surface salinity (SSS). Thirdly, we used passive tracers that behave in the model as any buoyancy-free particle or passive chemical from the rivers. Several experiments led to an identification of the river plume as the area where the tracer concentration exceeds 7 unit/m 3 , which best fits with the area identified by the criterion on SSS in the low runoff period (dry season).

The maps of the plume area show that south of 20°N, the RR plume is quickly joined by the plumes of three rivers south of the delta (Ma, Yen, Lam, hereafter the 'southern rivers'), creating a unique buoyant plume which extends southward along the coast most of the year. Therefore, we analyzed the variability of the resulting plume from the RR and the three southern rivers together. We first computed the plume surface area (in km 2 ) and found that the plume area follows the same variations as the river runoff, albeit with a time lag. Due to the interannual variability of the runoff, the plume area varies significantly between different years. Both the total discharge and the plume area undergo a strong variability at shorter time scales of a few days as well.

RR plume variability

The temporal variation and the regimes of the plume surface patterns.

The variation of daily spatial patterns of the plume is examined using KMA. The first cluster usually occurs from October to March, and is the most frequent cluster. In this cluster, the plume is very narrow and is mostly confined to the shallow area (depth < 20 m). The second cluster appears on and off throughout the year, with a slightly greater rate of occurrence in April and September which is during the seasonal transition of the monsoon. In this cluster, the river plume extends both further offshore and further southward compared to cluster 1. The third cluster happens primarily in June and July. In this cluster, the river plumes from the various rivers are disjoint. The fourth cluster exhibits the largest spatial coverage and occurs mostly in August and September. In this cluster, the plume extends the farthest offshore. Each cluster is associated with different forcing conditions (runoff, wind, current) that explain its spatial and temporal variability.

The thickness of the plume.

From September to March, the plume is narrow, elongated along the coast and attached to the bottom due to the winter monsoon downwelling winds. From May on, the monsoon changes to southwesterly wind, driving the plume northward and detaching it from the coast and from the bottom. In that case, the plume is surface advected with a thickness of 2-3 to 7-8 m. We observe that events where the plume detaches from the coast at 19°N coincide with the temperature at the coast colder than seaward. This supports the assumption that the southerly monsoon winds generate a coastal upwelling locally. Sometimes in May-August, the plume is thickening offshore to 12-15 m. That deepening may be due to strong winds events that increase the mixing offshore in the upwelling region or to the increase of river discharge.

In mid-August, we observed a strong plume deepening over the shelf (106-106.5°E). The deepening events appear to be linked with a seasonally recurrent anti-cyclonic eddy developing near 19°N. In our simulation, this eddy generally occurs in August, when the wind direction is from the southwest, which is upwelling-favorable, and the runoff is high. Then, this eddy disappears when the coastal southward current develops again and/or the wind is not upwellingfavorable anymore. The formation mechanism of this eddy is still unclear. In GOT_NORIV, some vorticity gradients can be depicted, suggesting an anticyclonic structure, but its amplitude is much weaker than in the reference simulation and it vanishes on a shorter time scale. The anticyclone seems strongly connected to the Lam river plume. Bottom topography may be another forcing factor since the eddy is located in the vicinity of an island.

Impacts of tides.

The influence of tides in shaping the far field plume is investigated by comparing the reference simulation GOT_REF with the simulation without tides (GOT_NOTIDE). On average, the plume in GOT_REF is 5% larger than in GOT_NOTIDE, with the largest differences found in the high discharge period. The KMA analysis on the GOT_NOTIDE shows that the 4 clusters are very close to the ones from GOT_REF, both in terms of spatial structures and temporal distribution.

However, tides impact the vertical structure of the plume: they enhance the vertical mixing hence the riverine water concentration in the bottom layer. They also lead to a larger export offshore of the riverine water at the surface close to the mouth during the high discharge season.

As this study focuses on the far field, we have not looked into details at the impact of tides on the water properties and dynamics at the mouth of the rivers, through the tidal-induced bottom friction. This should be done in future work.

Model response due to perturbations of the wind forcing

By comparing the daily ECMWF fields to the ASCAT data and ECMWF climatological mean, the uncertainty of the ECMWF wind in this area was estimated as 60% of the wind variability. The wind perturbations were then created as linear combinations of bivariate wind EOFs weighted using a Gaussian pseudorandom number generator with standard deviation equal to the uncertainty value estimated above. An ensemble of SYMPHONIE simulations was run with 50 members, perturbed by pairs of opposite perturbations so that the ensemble average of perturbations is zero.

The assessments are mostly based on the ensemble spread, which is computed as the standard deviation along the ensemble with respect to the ensemble mean.

The spread of the model variables.

The results show that the SSH spread is less than 2 cm most of the time. The mean spread of the SSH over the whole period decreases from the north to the south. It increases locally near the coast where the river plume is present. The peak of the spread occurss at the same time as the peaks of the wind stress spread, which is also the peak of the EOF1, and corresponds to a reversal of the monsoon over a few days..

The surface currents are highly sensitive to the wind perturbations, as expected. In particular, the high spread along the coast in V component suggests that the wind perturbation has an impact on the coastal current created by the river runoff. The coastal current can have opposite direction between two members at daily time scales, while the spread can exceed the current amplitude of the ensemble mean.

The temporal evolution of the spatial mean spread of SST shows the same time pattern as the SSH. However, unlike the SSH, the maximum spread of SST is highest near the south of Hainan island. We expect that the wind perturbations lead to changes in vertical mixing (which is already high due to tidal effect in that area) and then leading to the displacement of the water mass that has a high temperature gradient. Over the whole period, the mean spread is highest over the Red River coast, which can be related to the difference between temperature of the riverine water or the upwelling present in this area. The SSS, as expected, has high spread near the Red River delta.

The vertical section at 19°N is investigated to assess the vertical structure of the spread of temperature and salinity. The temperature spread is highest around 20-30m deep, where the vertical gradient is highest (base of the mixed layer; as in [START_REF] Andreu-Burillo | Assessing oceanmodel sensitivity to wind forcing uncertainties[END_REF][START_REF] Ayoub | Estimating uncertainties on a Gulf Stream mixed-layer heat budget from stochastic modeling[END_REF], a clear confirmation of the fact that the uncertainty of the wind affects vertical mixing. On salinity, the spread is highest near the surface in the vicinity of the sanity gradient due to the runoff.

The spread of the Red River plume.

At the surface, the spread of the plume area is highest in August, but it does not follow the same trend as the SSS. The maximum spread of the plume area happens at the same time as the maximum of the plume area, suggesting that our results might be sensitive to the uncertainty of the river runoff as well.

The clustering analysis was carried out per-member using as a reference the clusters of GOT_REF, as computed in Nguyen Duy et al. (2021). Results evidence some changes of clusters across the ensemble, but almost all the time, the pattern that is most likely to occur is still the reference cluster.

Finally, the spread of the plume thickness shows that the plume is not only changed vertically but also changed horizontally. In June, the spread is less than 1m. The peak of the spread happens from 09-13 July (> 3m), which is the same time as for the largest spread of wind stress, confirming the impact on vertical mixing. The plume can deepen inside eddies, as found near 19°N and as discussed in [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF]. The wind perturbation can alter the current, change the shape of the eddies and finally modify the plume thickness.

The spread of the plume for some specific days is also examined. The distribution over the ensemble of the vertical mixing coefficient evidences the increase of mixing in members with stronger wind. The mixing can in turn lead to a decrease of the tracer concentration at surface; if the concentration decreases below the threshold chosen to identify the plume, then the plume is not identified anymore, while a thicker mass of water is influenced by the river input. This suggests that using two conditions to identify the plume (one based on the surface concentration and the other one of the thickness of the water column impacted by the riverine waters) could be more robust to uncertainties in the wind forcing. Stronger winds also favor a spreading of the plume further offshore. However, disentangling the wind effect on vertical mixing and on advection is not straightforward. First, because the horizontal transport is also driven by the plume buoyancy, at least in the mid field (Horner-Devine et al., 2015). Second, because the ocean conditions at a given day result from the cumulative effect of forcing over a 'certain' period before (maybe several days). This should be taken into consideration when comparing the plume properties between different members at a given date, but as the daily variability of the wind is quite high in the GOT, such an approach is quite complex.

The ensemble has been assessed empirically using the SST and HFR data. In both cases, the ensemble is underdispersive, with more value fall in the first and last rank of the rank histogram. However, the error budget analysis shows that the difference between the ensemble spread and the ensemble-mean RMSE is compatible with the HFR observational error estimates.

Finally, to conclude, we estimate that within our experimental conditions the findings of [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF] are yet unchallenged and robust with respect to wind perturbations.

VI.2 Perspective

Model configurations and assessments

The assessment of the model representation of the plume dynamics evidences a lack of systematic high-resolution observations of the plume area, except for the limited area covered by the HF radar. The contribution of present altimetric measurements from the Jason satellite series and Sentinel-3 satellites should be investigated, but their spatial coverage (limited to the satellite nadir) may still be an issue. The future SWOT mission, with wide-swath altimetry (e.g. [START_REF] Fu | SWOT: The Surface Water and Ocean Topography Mission: Wide-Swath Altimetric Measurement of Water Elevation on Earth[END_REF] will represent a large amount of additional data which are expected as a precious contribution to both the delta and plume observation.

In this study, the river configuration is set up as realistic as possible, however, in the future, further improvements can be made.

The first improvement can be related to the runoff data. In this study, the total runoff in the Red River is totally the runoff data from the Son Tay hydrological station. However, there are other sources of the runoff, from Cau, Thuong, and Luc Nam river, that were not taken into account in this study because the discharge data was not available to us. The Web portal of Bac Giang (2016) indicates that the total runoff from 3 rivers is 7.52 billion m 3 /year, which is equal to around 238 m3/s while Vinh et al. (2014) estimates a different value of 133.5m 3 /s. These values are equivalent to around 4 -7% more of the total runoff in the RR. Also, the ratio of the runoff between the high discharge and low discharge period for the southern rivers is higher than RR. This can be explained by the fact that the RR runoff is controlled by dams. However, recently, some dams were also constructed in the southern rivers (e.g. Trung Son dam commissioned in 2017). This should be considered if the simulation period is extended in the future.

The second improvement can be the better distribution of the runoff over the RR mouth. In this study, the RR runoff is distributed throughout different distributaries using the annually averaged percentage of discharge made by Vinh et al. (2014). However, as indicated in their study, this percentage also has seasonal variation. For example, in Van Uc river, the discharge percentage for dry and rainy season is 14.1% and 14.6%, while for Cam, it is 9.9% and 9.4%, respectively. Therefore, another ensemble study could be useful to assess the impact of the uncertainty of the discharge rate and the distribution percentage on the plume variability. [START_REF] Warrick | Coastal river plumes: Collisions and coalescence[END_REF] suggested that more studies on interactions among different plumes are needed in general. As shown by [START_REF] Gong | Plume-to-plume interactions in the Pearl River Delta in winter[END_REF], a plume can have an impact on another plume, therefore the interaction between plumes is important for environmental management. With different passive tracers, the interaction between different rivers in the RR plume can be examined and can be a promising topic to explore.

RR plume variability

To identify the main regime of the surface river plume, K-means clustering method is used in this study. Though 4 clusters is enough to describe the main pattern, it is shown that the standard deviation of the SSS inside cluster 3 and cluster 4 is higher than inside cluster 1 and cluster 2, suggesting that the forcing conditions in clusters 3 and 4 are more unstable. A nested clustering inside these clusters can be done to further explore the plume pattern at the smaller scale.

The impact of tides on the shape of the plume has been explored, however, an in-depth analysis of the tidal processes, of the tides-river flow interactions and of their impact on the velocity, shear and stratification in the estuaries and at the mouths is still an open question. It will be the topic of a future study aimed at better understanding the shape and dynamics of the outflow in the near and mid fields. The dynamics at the mouth is also highly important for sediment dynamics close to the mouth area.

The dilution of river waters due to stirring by eddies or mixing could be explored further than what we started in this study. Actually, hourly outputs of currents from the GOT_REF and GOT_noriv simulations have been used by Cuong Manh Tran and Alexei Sentchev (LOG) to compute Lagrangian diagnostics (such as diffusivity) and compared them with diagnostics from HF radar currents. An article by Cuong Manh Tran is in preparation and should be submitted in March 2022 ('Assessment of relative dispersion in the Gulf of Tonkin using numerical modeling and HF radar observations'); I am co-author of this article.

The study of mixing due to instabilities at the front of the plume (eg Horner-Devine et al., 2015) or within the plume (such as those described by Ayouche et al. 2021 in the Gironde River plume for instance) was beyond the scope of this work. However, such submesoscale processes should be quantified in dedicated studies to better understand the dynamics at work in the Red River plume and to better anticipate the modelling and observational needs for an accurate and reliable monitoring and forecasting system of the plume. Similarly, the impact of surface waves (swell or sea state) on the mixing and their interactions with the current within the plume are other processes that would deserve further work.

Model response due to perturbations of the wind forcing

In this study, we did not take into account the spatial structure of ASCAT errors, because we do not know those errors in the GOT region. As a consequence, we limited ourselves in this study to using the comparison between ECMWF and ASCAT to get a very rough estimate of ECMWF error variances. Future studies should explore methods to extract more information from the satellite wind to generate the perturbations.

We also limited the source of stochastic perturbations to the surface wind only. Although several studies, cited in the introduction of this Chapter, agree on the wind being an essential, if not dominant, source of ocean uncertainties, we could have considered other sources. For instance, we became aware during this study that the runoff could be another major source of uncertainty impacting the RR plume, along with the wind. As another example, we could not fully clarify the apparent partial correlation between the spread of surface currents and the SSH spread -in order to better explore that link, atmospheric pressure would need to be perturbed in a consistent way with the wind, in particular since the GoT is a shelf sea and semi-enclosed sea, so it should be sensitive to water displacements as, e.g., a low-pressure center passes through. One way to do this comparison in a meaningful way would be to use multivariate ECMWF ensembles, perturb both the wind and the pressure, and perhaps add a wave model to better model the wind effects on currents. We did not have the time to explore those two directions, but it will need to be done.

Due to the constraint in computing resources, in this study, the ensemble has 50 members. This value is based on our experience in previous studies with the similar configurations(for example [START_REF] Vervatis | Assessment of a regional physical-biogeochemical stochastic ocean model. Part 1: Ensemble generation[END_REF]) with 50 members, Quattrocchi et al. (2014) with 40 members, [START_REF] Vervatis | Assessment of a regional physical-biogeochemical stochastic ocean model. Part 1: Ensemble generation[END_REF] with 40 members). Therefore, a sensitivity test on the impact of the ensemble size on the spread estimate should be done in the future.

In addition to the analysis of the ensemble spread, another approach to characterize the ensemble could be based on the clustering analysis. This would allow us to assess the pattern variability across the ensemble, for example, clustering across the members for a given day to assess the general regime of the model response.

We could also compute ensemble covariances in order to better characterize the spatial structures of the model uncertainties and to identify the relationships between the uncertainties on different variables (see for instance [START_REF] Echevin | Horizontal and vertical structure of the representer functions for sea surface measurements in a coastal circulation model[END_REF].

Finally, ensemble empirical assessment can be carried out to assess the realisticness of the ensemble.

sans marées (GOT_NOTIDE) est réalisée pour évaluer l'impact des marées sur les principaux modes de variabilité du panache fluvial.

2. Le panache de la rivière a été identifié à l'aide de différentes méthodes, ce qui nous a permis de sélectionner la plus appropriée compte tenu des objectifs. La méthode sélectionnée est basée sur des traceurs passifs qui sont injectés dans les rivières. En utilisant les outils décrits ci-dessus, combinés aux observations, nous avons obtenu les principaux résultats ci-dessous :

Evaluation du modèle

La simulation GOT_REF est évaluée en la comparant aux observations disponibles : données satellitaires (altimétrie pour la marée, OSTIA pour la température de surface de la mer, SMOS pour la salinité de surface de la mer), données in situ (VITEL, CFO, avec des profils de température et de salinité dans ces deux jeux de données), et radar HF (pour le courant de surface). Nous comparons également nos résultats avec l'OGCM (le modèle de circulation global utilisé aux frontières ouvertes) pour déterminer la valeur ajoutée de notre modèle régional à haute résolution par rapport à un modèle global. Globalement, malgré quelques biais, le modèle reproduit bien les conditions saisonnières, même dans la zone très peu profonde. La comparaison avec les données CTD près de la côte montre une surestimation de la salinité de surface, donc nous soupçonnons que les débits prescrits peuvent être sous-estimés.

Nous discutons également l'échelle de longueur horizontale typique (premier rayon de Rossby barocline) calculée à partir du champ de densité simulé et décrivons la circulation générale.

Identification du panache du Fleuve Rouge

Pour identifier le panache du FR, différentes méthodes sont envisagées. Tout d'abord, nous appliquons une méthode basée sur l'indice de stratification (SI), en supposant que l'eau du fleuve forme une couche qui flotte au-dessus de l'eau de l'océan. La deuxième méthode consiste à identifier le panache à partir de la salinité de surface de la mer (SSS). Dans la troisième, nous utilisons des traceurs passifs qui se comportent dans le modèle comme des particules sans flottabilité ou des produits chimiques passifs provenant des rivières. Plusieurs expériences ont permis d'identifier le panache fluvial comme la zone où la concentration du traceur dépasse 7 unités/m 3 , ce qui correspond le mieux à la zone identifiée par un critère en SSS dans la période de faible débit (saison sèche).

Les cartes du panache montrent qu'au sud de 20°N, le panache du FR est rapidement rejoint par les panaches de trois rivières situées au sud du delta (Ma, Yen, Lam, appelées par la suite "les rivières du sud"), créant un panache flottant unique qui s'étend vers le sud le long de la côte la plupart de l'année. Par conséquent, nous analysons la variabilité du panache résultant du FR et des trois rivières du sud ensemble. Nous calculons d'abord la surface du panache (en km 2 ) et trouvons que la surface du panache suit les mêmes variations que le débit, bien qu'avec un décalage dans le temps. En raison de la variabilité interannuelle du débit, la surface du panache varie significativement entre les différentes années. Le débit total et la zone du panache subissent également une forte variabilité à des échelles de temps plus courtes de quelques jours.

Variabilité du panache du Fleuve Rouge

Les variations temporelles et les régimes du panache en surface La variabilité saisonnière du panache est examinée à l'aide de l'analyse en K-means et de scènes quotidiennes du panache. Le premier cluster se produit habituellement d'octobre à mars et est le cluster le plus fréquent. Dans ce groupe, la forme du panache est très étroite et le panache est principalement confiné à la zone peu profonde (profondeur < 20 m). La faible différence entre les différentes formes que revêt le panache dans les différentes scènes associées au cluster 1 signifie que la condition de forçage est relativement stable. Ce cluster est associé à la période de débit le plus faible, sous un fort vent de nord-est (mousson d'hiver), qui favorise le downwelling. Le courant de surface côtier est orienté vers le sud et atteint sa plus grande intensité dans le cluster 1, ce qui est cohérent avec le vent. Toutes ces conditions favorisent le confinement des eaux faiblement salées du panache (SSS <30) sur la côte et leur extension vers le sud tout au long de la côte.

Le deuxième cluster apparaît par intermittence tout au long de l'année, avec un taux d'occurrence légèrement plus élevé en avril et en septembre (généralement avant l'occurrence du cluster 3 et après celle du cluster 4), ce qui correspond à la transition saisonnière de la mousson. Ceci suggère que le cluster 2 représente un régime de transition pour le panache. Dans ce cluster, le panache fluvial s'étend à la fois plus au large et plus au sud que dans le cluster 1. Dans le troisième cluster, qui se produit principalement en juin et juillet, le panache a une forme différente de celle des clusters 1 et 2; les panaches des différentes rivières sont disjoints. Le quatrième cluster se produit principalement en août et en septembre, c'est-à-dire au plus fort de la saison de fort débit. Le panache présente alors la plus grande couverture spatiale et s'étend le plus loin au large.

Chaque cluster est associé à différentes conditions de forçage (débit des rivières, vent, courant) qui expliquent leur variabilité spatiale et temporelle.

L'épaisseur du panache. Impacts des marées.

L'influence des marées sur le panache dans le far-field est étudiée en comparant la simulation de référence GOT_REF avec celle sans marées (GOT_NOTIDE). En moyenne, le panache dans GOT_REF est 5% plus grand que dans GOT_NOTIDE, avec les plus grandes différences trouvées dans la période de haut débit. L'analyse KMA sur la simulation GOT_NOTIDE montre que les 4 clusters sont très proches de ceux de GOT_REF, tant en termes de structures spatiales que de distribution temporelle.

Cependant, les marées ont un impact sur la structure verticale du panache : elles augmentent le mélange vertical et donc la concentration d'eau d'origine fluviale dans la couche inférieure. Elles conduisent également à une plus grande exportation au large de l'eau fluviale à la surface près de l'embouchure pendant la saison de haut débit.

Comme cette étude se concentre sur le far-field, nous n'avons pas examiné en détail l'impact des marées sur les propriétés hydrologiques et sur la dynamique à l'embouchure des rivières, par le biais de la friction du fond induite par les marées. Ceci devrait être fait dans un travail futur.

Réponse du modèle due aux perturbations du forçage du vent

En comparant les champs quotidiens ECMWF aux données ASCAT et aux moyennes climatologiques de l'ECMWF, l'incertitude du vent ECMWF dans cette zone est estimée à 60% de la variabilité du vent. Les perturbations du vent sont ensuite créées comme des combinaisons linéaires d'EOFs bivariés de vent, pondérées à l'aide d'un générateur de nombres pseudoaléatoires gaussien dont l'écart-type est égal à la valeur de l'incertitude estimée ci-dessus. Un ensemble de simulations SYMPHONIE est réalisé avec 50 membres, perturbés par des paires de perturbations opposées de sorte que la moyenne des perturbations de l'ensemble soit nulle.

L'analyse des ensembles en termes de sensibilité a consisté principalement à calculer l'écarttype de l'ensemble par rapport à la moyenne de l'ensemble (la "dispersion") des variables océaniques et les distributions statistiques des échantillons de l'ensemble.

La dispersion de différentes variables du modèle

Les résultats montrent que la dispersion de la hauteur de la surface de la mer (SSH) est inférieure à 2 cm la plupart du temps. La dispersion moyenne de la SSH sur l'ensemble de la période diminue vers le sud, mais augmente aussi localement près de la côte vietnamienne dans la zone du panache. Le pic de la dispersion en SSH se produit en même temps que le pic principal de la dispersion de la tension de vent. Il s'agit de la periode où l'EOF1 atteint sa plus forte amplitude et correspond à une renverse de la mousson sur quelques jours.

Les courants de surface sont très sensibles aux perturbations du vent, comme on pouvait s'y attendre. En particulier, la forte dispersion le long de la côte de la composante V suggère que la perturbation du vent a un impact sur le courant côtier créé par les eaux fluviales. Le courant côtier peut avoir une direction opposée entre deux membres à des échelles de temps quotidiennes, tandis que la dispersion peut dépasser l'amplitude du courant de la moyenne d'ensemble.

La moyenne spatiale de la dispersion de la température de surface de la mer (SST) montre les mêmes variations temporelles que la SSH. Cependant, contrairement à la SSH, le maximum de dispersion de la SST se produit près du sud de l'île de Hainan. Nous pensons que les perturbations du vent entraînent une modification du mélange (qui est déjà élevé en raison de l'effet de marée dans cette zone) et conduit ensuite au déplacement de la masse d'eau qui présente un fort gradient de température. Sur l'ensemble de la période, l'écart moyen est plus élevé sur la côte Vietnamienne, ce qui peut être lié d'une part aux déplacements du gradient de température entre les eaux fluviales et les eaux océaniques côtières et d'autre part à un upwelling présent dans cette zone. La SSS, comme prévu, présente une forte dispersion près du delta du Fleuve Rouge.

Nous évaluons la structure verticale de la dispersion de la température et de la salinité à partir de sections verticales à 19°N. La dispersion de la température est la plus élevée autour de 20 -30m de profondeur, là où le gradient vertical est le plus élevé à la base de la couche de mélange (comme dans [START_REF] Andreu-Burillo | Assessing oceanmodel sensitivity to wind forcing uncertainties[END_REF][START_REF] Ayoub | Estimating uncertainties on a Gulf Stream mixed-layer heat budget from stochastic modeling[END_REF]. Ceci est une confirmation claire du fait que l'incertitude du vent affecte le mélange vertical. Sur la salinité, l'écart est le plus élevé près de la surface, à proximité du gradient de salinité dû aux eaux des rivières.

Dispersion des caractéristiques du panache.

En surface, la dispersion de la surface de la zone du panache est la plus élevée en août, mais elle ne suit pas la même tendance que la dispersion en SSS. La dispersion est maximale quand la surface du panache est elle-même maximale, ce qui suggère une sensibilité significative de nos résultats aux incertitudes sur les débits également.

Une analyse en clusters est effectuée sur chaque membre, en prenant pour référence les clusters de GOT_REF, calculés et présentés dans [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF]. Les résultats mettent en évidence certains changements de clusters dans l'ensemble, mais presque tout le temps, le cluster le plus susceptible de se produire reste celui comme dans le run de référence.

Enfin, la dispersion de l'épaisseur du panache montre que le panache n'est pas seulement impacté verticalement mais aussi horizontalement par les perturbations en vent. En juin, la dispersion est inférieure à 1m. Le pic de la dispersion se produit du 09 au 13 juillet (> 3m), ce qui correspond à la plus grande dispersion de la tension de vent, confirmant l'impact des perturbations en vent sur le mélange vertical. Le panache peut s'approfondir au centre de tourbillons, comme dans ceux simulés près de 19°N et discutés dans [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF].

Les perturbations du vent peuvent modifier le courant, changer la forme du tourbillon et enfin modifier l'épaisseur du panache.

La dispersion d'ensemble du panache est également examinée. pour certains jours spécifiques. La distribution sur l'ensemble du coefficient de mélange vertical met en évidence l'augmentation du mélange dans les membres avec un vent plus fort. Le mélange peut à son tour conduire à une diminution de la concentration du traceur à la surface ; si la concentration diminue en dessous du seuil choisi pour identifier le panache, alors le panache n'est plus considéré comme présent, alors qu'une masse d'eau plus épaisse est influencée par l'apport fluvial. Ceci suggère que l'utilisation de deux conditions pour identifier le panache (l'une basée sur la concentration en surface et l'autre sur l'épaisseur de la colonne d'eau impactée par les eaux fluviales) pourrait être plus robuste aux incertitudes du forçage en vent.

Des vents plus forts favorisent également une extension du panache plus au large. Cependant, il n'est pas simple de démêler l'effet du vent sur le mélange vertical et sur l'advection. Premièrement, parce que le transport horizontal dépend aussi de la flottabilité du panache, au moins dans le mid-field (Horner-Devine et al., 2015). Deuxièmement, parce que les conditions océaniques à un jour donné résultent de l'effet cumulatif du forçage sur une " certaine " période antérieure (peut-être plusieurs jours). Cela devrait être pris en compte lors de la comparaison des propriétés du panache entre différents membres à une date donnée, mais comme la variabilité quotidienne du vent est assez élevée dans le GOT, une telle approche est assez complexe à mettre en oeuvre.

Finalement, en conclusion, nous estimons que dans nos conditions expérimentales, les résultats de [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF] ne sont pas remis en cause et sont robustes par rapport aux incertitudes sur le vent.

Figure I. 1 Figure I. 3 Figure I. 4 25 Figure I. 7 Figure II. 1 Figure II. 4

 13425714 Figure I.1 Concept of a prototypical plume with all the dynamical regions, from Horner-Devine et al. (2015). ........................................................................................................................................ 19 Figure I.2 (A) Location of the Gulf of Tonkin; the blue dots indicate the mouths of the Red River delta. Contour lines indicate the bathymetry of the area. (B) Same as A, but zoom to the blue box. The star indicates the location of Son Tay hydrological station at the apex of the delta, while the square boxes show the river mouths of the Red River (red square) and of other rivers to the south (light blue) and to the north (yellow). ............................................ 21 Figure I.3 Monthly averages of wind direction (arrow) and velocity (m/s, in colors) in GOT, calculated from ECMWF analysis in the period from 2009 to 2018. .............................. 22 Figure I.4 Daily averaged precipitation in different months, calculated from ECMWF analysis in the period from 2009 to 2018. Units are mm/day. .................................................... 22 Figure I.5 Co-tidal charts for the constituents of (a) K1, (b) O1, (c) M2, and (d) S2, adapted from Minh et al. (2014). (Solid line: phase-lag in degree, dashed line: amplitude in cm). .................................................................................................................................................................................. 24 Figure I.6 Schematic representation of the Red River delta and its main tributaries (from Vinh et al., 2014). .............................................................................................................................................. 25 Figure I.7 Signature of the plume from the Red River, as depicted in Landsat-8 data on 25 August 2014 using band 5, 3 and 6. The star indicates the location of Son Tay hydrological station. Landsat-8 image courtesy of the U.S. Geological Survey. Landsat Scene Identifier: LC81260462014268LGN01. ..................................................................................... 27

Figure I. 1

 1 Figure I.1 Concept of a prototypical plume with all the dynamical regions, from Horner-Devine et al. (2015).

  divides the prototypical plume into two parts: the near field where the outflow is critical and the advection and shear mixing dominate the dynamics; and the far field where the earth's rotation, wind, and background flow dominate the evolution of the plume. Horner-Devine et al. (2009), in their conceptual model, and Horner-Devine et al. (2015), in their review, describe the Columbia river in terms of four water masses: source water, the near field, the mid-field, and the far field (Figure I.1).

  Gulf of Tonkin (hereafter GOT) (16.9°N -21.9°N, 105.6°E -110°E) is a small shelf sea located east of Vietnam and south of China (Figure I.2A). It is a shallow area with depth less than 100m (Figure I.2B). It connects with the open ocean through the Hainan strait (with depth approximately 40 -50m) and south of Hainan (with depth approximately 80 -100m).

Figure I. 2 (

 2 Figure I.2 (A) Location of the Gulf of Tonkin; the blue dots indicate the mouths of the Red River delta. Contour lines indicate the bathymetry of the area. (B) Same as A, but zoom to the blue box. The star indicates the location of Son Tay hydrological station at the apex of the delta, while the square boxes show the river mouths of the Red River (red square) and of other rivers to the south (light blue) and to the north (yellow).

Figure I. 3

 3 Figure I.3 Monthly averages of wind direction (arrow) and velocity (m/s, in colors) in GOT, calculated from ECMWF analysis in the period from 2009 to 2018.

Figure I. 4

 4 Figure I.4 Daily averaged precipitation in different months, calculated from ECMWF analysis in the period from 2009 to 2018. Units are mm/day.

Figure I. 5

 5 Figure I.5 Co-tidal charts for the constituents of (a) K1, (b) O1, (c) M2, and (d) S2, adapted from Minh et al. (2014). (Solid line: phase-lag in degree, dashed line: amplitude in cm).

Figure I. 7

 7 Figure I.7 shows a Landsat-8 image, with the presence of the large spread of a water region with different optical properties near the coast created by the runoff from the RR. It is an example of the signature of the river plume.

Figure I. 7

 7 Figure I.7 Signature of the plume from the Red River, as depicted in Landsat-8 data on 25 August 2014 using band 5, 3 and 6. The star indicates the location of Son Tay hydrological station. Landsat-8 image courtesy of the U.S. Geological Survey. Landsat Scene Identifier: LC81260462014268LGN01.

Figure II. 1

 1 Figure II.1 The Arakawa staggered C-grid[START_REF] Piton | Du Fleuve Rouge au Golfe du Tonkin: dynamique et transport sédimentaire le long du continuum estuaire-zone côtière[END_REF]. The i, j, k indices correspond to x, y, z direction, respectively. φ represents the scalar fields. Turbulent fields are defined at the same positions as the vertical component of velocity.

Figure II. 2 (

 2 Figure II.2 (A): The grid size of GOT_REF (m). (B): Sigma vertical layers coordinate along the red transects in (A).

Figure II. 3

 3 Figure II.3 The river input points of different rivers: Red River (red), southern rivers (cyan), northern rivers (yellow). Red lines illustrate the channels used for the RR system. The river runoff flows from the most upstream point of the channel. Note: the length of the red lines in this figure does not represent the real lengths of the channels.

Figure II. 4

 4 Figure II.4 Daily Red River discharge in the model configuration, from National Hydro-Meteorological Service data (red curve), monthly climatological discharge of the Southern Rivers (black curve), monthly climatological discharge of the Northern Rivers (green curve). Units are m 3 /s.

Figure II. 5

 5 Figure II.5 The scheme and configuration of a channel in the model.

Figure II. 6

 6 Figure II.6 illustrates the salinity profile in the Ninh Co channel, which has the lowest percentage of runoff over the 7 channels, on 03 March 2015. In the simulation with tide, at 0h, the salinity is present at around 35km from the river mouth (Figure II.6A), while at 12h, the propagation is 10km less (Figure II.6B), with the contour line of 0.0001.In the simulation without tide, the salinity intrusion seems more stable within 1 day. The salinity in both time frames shows strong stratification at around 25km from the river mouth. It also shows the changes of SSH over the day in GOT_REF, while in GOT_NOTIDE the change in SSH is much less significant. The important thing is that in both cases, with and without tide, the salinity at the river mouth is not zero and the implementation of the channel can help create better salinity mixing at the river mouth.

Figure II. 6

 6 Figure II.6 Salinity sections from GOT_REF in Ninh Co on 03 March 2015, at 0h (A), 12h (B). (C, D): same as (A, B) but in GOT_NOTIDE

II. 2

 2 .3.1 Light penetration in the first layers of the ocean In order to model correctly the ocean dynamics, using a suitable heat flux parameterization is very important. The change of heat flux can lead to change of temperature in the ocean and then change the mixing. The change of mixing can then impact the variations of the river plume.

Figure III. 2 (

 2 in Chapter 3), shows a tongue of cold water near the Hainan island, which is consistent with such hypothesis (see discussion in III.2.1.2).

Figure II. 7

 7 Figure II.7 An example of monthly averaged KD490, downloaded from Ifremer, in (A): March 2010, and in (B): September 2010. Lower values means that the ocean is clearer. Units in m -1 .

Figure II. 8

 8 Figure II.8 L2 (m) calculated from KD490 (in color). The isolines show the bathymetry.

  Run the model ○ Set up the environments (since the model has been run on different machines: HILO at USTH, NUWA at OMP, CALMIP at University of Toulouse). ○ Update the model source in case of using new functions and/or change the parameters that are hard-coded in the source. ○ Compile the model. ○ Update the model parameters depending on the purposes. ○ Run the model ○ Repeat the whole process with each simulation. • Assess the model with observational data ○ Gridded data: interpolate the model results to the observational data grid if the spatial resolution of the data is lower than the model result. ○ Data on points: find the location of the model grid points which are nearest to the location of the data; usually the model-data comparison is made for several neighboring model grid points in order to account for representativity errors due for instance to errors in bathymetry. ○ Do the calculations Table

Figure

  Figure II.9 (A) Model domain with bathymetry (m, in color). Blue boxes show the location of river mouths of Red River. Cyan and yellow boxes show the location of other river mouths in the south and north of RR, respectively. XUAN and DHOI are the locations of the HF radars. (B): Location of VITEL CTD stations. (C) Location of CFO CTD stations (Red dots). (D) Tidal data. Red dots with black numbers: Location of the tide gauges from which historical data of tidal constituents are available. Station 2 and 3 stay nearly at the same location but measurements were made over different periods of time. Black dots: Location of the tracks from the T/P and Jason series altimetric missions from which tidal constituents are estimated.

Figure II. 10

 10 Figure II.10 Process of clustering the river plume in GOT.

III. 2 . 1

 21 Figure III.1 shows the comparison between GOT_REF and SMOS sea surface salinity (SSS) data (please refer to Chapter 2, section II.3.4 for more information about SMOS). In the winter, from December to January, the space-time average of SMOS SSS is 32.78. OGCM (the global simulation that we use as boundary condition) and GOT_REF show mean values of 32.90 and 32.95, respectively: the mean bias between the data and models is therefore small. In all cases, there is a low salinity region near the Hainan strait, which is probably created by the runoff of the Pearl River in the north (outside the GOT_REF domain). However, in SMOS, this low salinity water extends to the south (Figure III.1A), while in OGCM and GOT_REF it extends to the north along the coast and is connected with the other river plumes (FigureIII.1C, E). The pathway shown in OGCM and GOT_REF is more consistent with the cyclonic current in the winter described in previous studies(Ding et al., 2013; Gao et al., 2014; Rogowski et al., 2019; Piton et al., 2021). It is also possible that the SSS is too high in the OGCM which is used for initial and open-boundary conditions, leading to the overestimated SSS with respect to SMOS data in the center of the basin.

Figure III. 1 (

 1 Figure III.1 (A): Sea surface salinity from SMOS satellite, averaged from December to January, from 2011-2016. (B): Same as (A) but averaged from June to August. (C, D): Same as (A, B) but for OGCM. (E, F): Same as (A, B) but for GOT_REF. (G, H): mean error of SMOS, calculated from the same period as (A, B).

Figure III. 2 (

 2 Figure III.2 (A): OSTIA sea surface temperature, averaged from December to January in the period of 2011-2016. (B): Same as (A) but for OGCM. (C): Same as (A) but for GOT_REF. (D, E, F): Same as (A, B, C) but calculated from June to August. Unit in °C.

Figure

  Figure III.3 (A): Daily mean SLA computed from December to February in the period 2011-2016. (B): same as (A) but for OGCM. (C): same as (A) but for GOT_REF. (D, E, F): same as (A, B, C) but computed from June to August.

Figure III. 4 (

 4 Figure III.4 (A): Vertical salinity sections from CFO CTD data (left) and model result (right) on 05 September 2015. (B): same but for temperature. (C, D), (E, F), (G,H): same as (A, B) but for 15 January 2016, 14 April 2016, and 25 July 2016, respectively.

Figure

  Figure III.5 (A): Bottom temperature (°C) of GOT_REF, spatially averaged where the depth is less or equal to 30m (blue line, left axis); difference between GOT_23 and GOT_REF (red line, right axis). (B): same as (A) but the averages are computed over the region deeper than 30m.

Figure III. 6

 6 Figure III.6 Surface temperature (°C) of GOT_REF over the whole area (blue line, left axis); difference between GOT_23 and GOT_REF (red line, right axis).

  February and from April to September over the period 2011-2016. In GOT_23, the MLD ranges from 2m near the coast to around 8m in the deep region in summer (Figure III.7A), while it reaches 30 to 50m in the deep region in winter (not shown).

Figure

  Figure III.7 (A): Mixed layer depth of GOT_23, calculated from the mean density profile from April to September in the period 2011-2016. (B): same as (A) but for GOT_REF. (C): the difference between (A) and (B). Units in m.

Figure

  Figure III.8 (A): Temperature profile from VITEL CTD data. Temperature profile from GOT_23 (B) and from GOT_REF (C). Unit in °C.

  Figure III.9 shows the mean of the different heat flux components over the whole domain in different months, calculated from 2011-2016.

Figure

  Figure III.9..(A) Total heat fluxes in W/m 2 from GOT_23 and GOT_REF, calculated from 2011-2016, averaged over the whole domain. Black line shows the solar flux (solar flux is the same for GOT_REF and GOT_23). (B): Latent, sensible and longwave heat flux from GOT_23 (left column) and GOT_REF (right column). Positive (resp negative) values mean a warming (resp cooling) of the ocean..

Figure

  Figure III.10 RMS of the difference between daily salinity fields from GOT_REF and from GOT_23 at different depths computed over the whole period 2011-2016.

FIGURE 1 |

 1 FIGURE 1 | (A) Location of the Gulf of Tonkin. Blue dots show the location of the river mouths of the Red River system. Contour lines indicate the bathymetry of the area. (B) Model domain with bathymetry. The star indicates the location of Son Tay hydrological station at the apex of the delta. Blue boxes show the most upstream points where the runoffs flow to the model. Cyan and yellow boxes show the location of other river mouths in the south and north of RR, respectively. (C) Daily Red River discharge in the model configuration, from National Hydro-Meteorological Service data (red curve, left y-axis). Blue dashed line: Percentage of Red River runoff compared to total runoff (sum of runoff from Red River, other rivers in Vietnam, rivers in China). Blue line: percentage of all rivers in Vietnam (including Red and other rivers in the GOT) compared to total runoff. Blue lines are referred to the right y-axis.

FIGURE 2 |

 2 FIGURE 2 | (A) Grid cell size (m). Red squares show the location of HFR antennas. (1, 2, 3): Sections used to calculate the tracer transport. VITEL (B) CTD stations location. Red dots indicate the station locations, while blue points indicate the river mouths' locations.

(

  NHMS) of Vietnam at Son Tay hydrological station which is located at the apex of the RR delta. The discharge is distributed across the seven distributaries based on the results ofVinh et al. (2014) (Bach Dang (7%), Cam (13.2%), Van Uc (14.5%), Thai Binh (6.4%), Ba Lat (30.3%), Ninh Co (5.6%), Day (23.0%)). Furthermore, each river mouth is connected to a channel in the model to allow the water entering the coastal ocean with realistic salt and temperature properties and realistic stratification. The length of the channels (35-45 km, depending on the channel) is chosen to exceed the saltwater intrusion, which is approximately 30 km from the mouth[START_REF] Nguyen Thi Hien | Upstream effects on salinity dynamics in the Red River Delta[END_REF]. The results from GOT_REF confirm that the salty water never reaches the end of the channel, even in the low discharge period. The river runoff is converted into a vertically sheared current at the most upstream point of the channel (APPENDIX B); there the salinity of the river flow is and the temperature varies seasonally from 17 • C (in February) to 29 • C (in August).

FIGURE 3 |

 3 FIGURE 3 | Amplitude of the tidal constituent (in colors, cm) from GOT_REF. Each figure shows the complex error (mm), calculated from the amplitude and phase, between model result and observational data for 4 main tidal constituents (O1, K1, M2, S2). Black circles show errors compared to altimetry data, while red circles indicate errors compared to the tidal gauges data.

FIGURE 4 |FIGURE 5 |

 45 FIGURE 4 | (A,B) The first baroclinic Rossby radius of deformation (km), over the period 2011 -2016, computed from the temperature and salinity profiles averaged from April to September (A) and October to March (B). Contour lines indicate the bathymetry. Star shows the location of Bach Long Vi island. Note: Different color bar scales.

FIGURE 6 |FIGURE 7 |

 67 FIGURE 6 | Surface current from HFR data (A,C,E,G,I) and model simulations (B,D,F,H,J) for different months of 2015 (Jan, Apr, May, Jul, Oct) (cm/s). Temporal correlation of u and v (K,L) computed from daily currents. Note: Different color bar scales.

  use SSS = 29 to detect the plume, whileBurla et al. (2010) use the value of SSS = 28 and in[START_REF] Maccready | A model study of tide-and wind-induced mixing in the Columbia River Estuary and plume[END_REF] the plume is identified using SSS < 31.

FIGURE 9 |

 9 FIGURE 9 | Frequency of occurrence of the plume within cluster 1 (A), cluster 2 (B), cluster 3 (C), cluster 4 (D). For instance, for cluster 1: at a given point, a frequency of 25 means that over the whole period when cluster 1 is present, the plume is present at this point 25% of the time. (E) Temporal distribution of each cluster (in red) and the runoff from Red and Southern rivers (in black).

  FIGURE 10 | (A-D) Wind conditions (m/s) corresponding to the 4 clusters. (E-H) and (I-L) Same as (A-D) but for surface current (cm/s) and SSS. (M-P) SSH anomaly with respect to the mean over the basin for each cluster (cm).

  FIGURE 11 | (A) Plume thickness (m) at 19 • N section from 2013 to 2014. Red points indicate the area where the plume is present over the whole depth. (B) Temperature ( • C) at 5 m at the same section as (A). (C,D) Tracer concentration over time at two coastal points. The location of the section and the points are shown in panel (E).

  FIGURE 12 | (A-C) Surface tracer concentration on 01 June 2012 (A), 07 June 2012 (B), and 16 June 2012 (C). (D,E) Relative surface vorticity (10 -6 s -1 ) (contour fill, left color bar) and plume thickness (m) (line, right color bar) on 11 and 15 August 2014. The black box and red circle indicate the location of two islands.

FIGURE 14 |

 14 FIGURE 13 | (A-D) Maps of mean surface tracer concentration correspond for 4 clusters. (E-H) Difference of surface tracer concentration between GOT_REF and GOT_NOTIDE correspond for 4 clusters. (I-L) Same as (A-D) but at the bottom. (M-P) Same as (E-H) but at the bottom.

Figure

  Figure S2: (A, B, C, D): Average salinity for GOT_REF for each cluster. Black contour indicates the 30 isohalines. (E, F, G, H): Same as (A, B, C, D) but for GOT_noriv. (I, J, K, L): Difference between GOT_REF and GOT_NORIV. Red and black contour lines indicate the difference in SSS of 0.5 and 1, respectively. Black dots show the area where the SSS difference is less than 0.1.

Figure S3 :

 S3 Figure S3: Plume thickness at 19 N (A) and at 20N (B) from 2011 -2016. Red points indicate the area where the plume is present over the whole depth.

Figure S4 :

 S4 Figure S4: (A, B, C): Relative surface vorticity (10 -6 s -1 ) (contour fill, right color bar) and surface current direction (arrow) in GOT_REF on 15, 20, and 22 August 2014, respectively. (D, E, F): same as (A, B, C) but for GOT_noriv.

Figure

  Figure IV.16 (A) Surface tracer concentration of GOT_REF (sponge layer = 30 grid points) on 15 August 2011. White dashed line shows the limit of the sponge layer. (B): Surface concentration at the same time but for a GOT_S10 (sponge layer = 10 points).

Figure

  Figure.IV.17 (A): Black line: depth-averaged current of section 1 (positive means northward). Blue line: Difference of depth-averaged salinity between GOT_REF and OGCM at blue point. (B): same as (A) but for section 2 and green point. (C): same as (A) but for section (3) and red point. The locations of section 1, 2, 3 and blue, green, red points are shown in (D).The salinity is taken at the circle for GOT_REF and at the square for OGCM. If the black line is positive and the colored line is negative, the unreal freshwater will be introduced to the model.

  us summarize what has been done so far. Using numerical results combined with observational data, the variability of the RR plume has been examined by several authors. Gao et al. (2013) and Rogowski et al. (2019) attempted to explain the plume variability in the summer and winter. The variability of the RR and southern rivers plume is further examined by Nguyen-Duy et al. (2021) using a high resolution, multiple sources cross validated bathymetry, pluri-annual 3D model. As shown in Nguyen-Duy et al. (

Figure V. 1

 1 Figure V.1 displays the tracer concentration, SSS, and wind conditions, corresponding to cluster 3 and cluster 4, but averaged from June to August 2015 instead of over the whole period(2011)(2012)(2013)(2014)(2015)(2016) as presented in Chapter 4. Though averaged over a much shorter period with respect to the published article, the figure still shows the characteristic behavior of the plume under different wind conditions. Under strong upwelling favorable wind, the plumes from the different rivers are disconnected, and advected more to the north (Figure V.1A), while they form a single plume when the upwelling wind is less intense (Figure V.1D). Secondly, regarding the plume thickness, in winter, when the wind is downwelling favorable, the plume sticks to the coast and reaches the bottom. In summer, due to monsoon, the wind is upwelling favorable and even if the runoff is three times higher than in winter, the plume is thinner and advected offshore.

Figure V. 1 (

 1 Figure V.1 (A): Surface tracer concentration (arbitrary unit/m 3 ), averaged from daily concentration, from June to August 2015 for the days when the plume is classified into cluster 3. (B): Same as (A) but for SSS. (C): Same as (A) but for 10m wind from ECMWF (in m/s). (D, E, F): same as(A, B, C) but averaged for the days when the plume is classified into cluster 4. The contour lines in (A, D) show the area where the surface concentration is equal or higher than 7 unit/m 3 , which is the threshold to be identified as the river plume.

Figure V. 2 (

 2 Figure V.2 (A) The spatial component of EOF 1. The colors and arrows present the surface wind velocity and direction, respectively. (B) The temporal component of EOF 1.(C, D), (E, F), (G, H), (I, J): same but for EOF 2, 3, 4, and 5, respectively. (K): Fraction of variability explained by the first 12 EOFs. The upper number is the ratio of the EOF variance to the total variance, while the lower number in the bracket is the cumulative variance.

Figure V. 3

 3 Figure V.3 (A): Standard deviation of the bias between ASCAT and ECMWF, for U component (ms-1). (B): Standard deviation of daily mean ECMWF wind for U component (ms-1). (C): ratio between (A) and (B). (D, E, F): same as (A, B, C) but for V component.

Figure

  Figure V.4A shows the temporal standard deviation of the bias for U between the daily ECMWF in 2015 and ECMWF_MEAN. The standard deviation is lower near the RR coast (1.7 ms -1 ) than at the south (3 ms -1 ). The ratio between the standard deviation of the bias and standard deviation of the wind is shown in Figure V.4C. It is lower in the south and the center of the gulf (~0.9) and higher in the north (~ 1 to 1.05).

Figure V. 4

 4 Figure V.4 (A): Standard deviation of the bias between the climatological ECMWF mean and ECMWF, for U component. (B): Standard deviation of daily mean ECMWF for U component. (C): ratio between (A) and (B). (D, E, F): same as (A, B, C) but for V component. Units for (A), (B), (D), (E): ms-1.

Figure

  Figure V.5 C, D shows the mean of ensemble spread of wind stress over the period, for U and V components, respectively. While the spread of the U component is higher in the south of the gulf, while the spread of V component is higher in the north. This is consistent with the standard deviation of the U and V components of ECMWF (Figure V.4 B, E). V.3.1.1.2 SSH

Figure

  Figure V.6 (A): Temporal evolution of daily spatial mean of the ensemble spread of sea surface height (SSH) for ENS_REF and ENS_0.3. (B): mean of daily SSH ensemble spread of ENS_REF over the whole period, excluding the period from 01 -15 July. (C): SSH spread of ENS_REF on 10 July 2015. Units are cm.

Figure V .

 V Figure V.6A shows the temporal evolution of the spatial mean of the ensemble spread of SSH for ENS_REF and ENS_0.3. Most of the time, the SSH spread < 2cm (for ENS_REF) and < 1 cm (for ENS_0.3). Overall, the mean spread of ENS_REF and ENS_0.3 are 1.48 and 0.71 cm,

  Figure V.7 shows the temporal evolution of the spatial mean of the ensemble spread of the surface current for ENS_REF. Most of the time, the spread of U and V components is similar, around ~ 4 cm s -1 , except for the period of 5 -11 July (mostly the meridional currents), which is the same time when the spread of wind and SSH increases. Some other higher spread events are 08 August and 21 August, which correspond also to high SSH spread.

Figure V. 7

 7 Figure V.7 Temporal evolution of daily spatial mean of the ensemble spread of surface current for ENS_REF (cm s -1 ).

  Figure V.8 C, D), which is even larger than the ensemble mean (5 cm s -1 for U and 10 cm s -1 for V component -Figure V.8 E, F).

Figure V. 8

 8 Figure V.8 Ensemble spread of surface current on 09 July 2015, for U component (A) and V component (B). (C, D): same as (A, B) but zoom on the red box. (E, F): same as (C, D) but for ensemble mean. Units are cm s-1.

Figure

  Figure V.9 (A): Temporal evolution of daily spatial mean of the ensemble spread of sea surface temperature (SST) for ENS_REF and ENS_0.3. (B): SST spread of ENS_REF on 11 July 2015. (C): SST of GOT_REF on the same day as (B). Units are °C.

Figure

  Figure V.10 (A): SST spread of ENS_REF on 11 July 2015, but zoom to RR coast. (B): SST of GOT_REF on the same day. Units are °C.

Figure V .

 V Figure V.11A shows the temporal variation of the SST spread at 20°N. The high spread tends to occur where the horizontal gradient of the SST is high (Figure V.11B). This gradient can be created by the upwelling in June and early July, as suggested by Nguyen-Duy et al. (2021) and can be seen here in the signature of the temperature at the surface, 5 m and 10 m of GOT_REF (Figure V.11 B, C, D). High spread can occur also when the wind and surface current spreads are large (July 8 -12). The largest spread occurs near the coast on 23 August, but there is no signature of the upwelling in the SST. Instead, it coincides with a local maximum of spread in SSS (not shown). As said above, it can be due to the difference between the temperature of the river and ocean water (Figure V.11 D).

Figure V. 11

 11 Figure V.11 (A): Hovmuler diagram of the SST spread at 20°N as a function of longitude and time. (B, C, D): GOT_REF temperature at the surface, 5m, and 10m, respectively. Units are °C.

Figure

  Figure V.12 (A): Temporal evolution of daily spatial mean of the ensemble spread of sea surface salinity (SSS) for ENS_REF and ENS_0.3. (B): Same as (A) but for surface tracer concentration of ENS_REF. Time-mean of daily SSS spread of ENS_REF over the whole period (C) and zoom to the RR coast (D).

Figure V .

 V Figure V.12Cshows the spatial distribution of the ensemble spread on 11 August when the spread is the highest. It shows that the spread is less than 1 in almost all areas, except near the coast: where the RR plume is present, it can reach 4. The SSS spread over the whole period shows the same pattern, with a high spread near the RR coast (Figure V.12C). Near the southern boundary, the SSS spread is less than 0.1. Values higher or equal to 0.5 are only reached near the coast near the river mouths.

  Figure V.13A. The maximum spread reaches 1.5°C at around 30m depth. The same section of GOT_REF shows that the maximum spread occursat the maximum gradient depth (base of the mixed layer) (Figure V.12B, C) consistently with the findings of Andreu-Burillo et al. (2002)

Figure

  Figure V.13 (A): Zonal ensemble spread of temperature at 19°N on 11 August 2015 (°C). (B): Vertical section on the same day off GOT_REF (°C). (C): vertical gradient of (B) (°C/m). (D, E, F): same as (A, B, C) but for salinity.

Figure

  Figure V.14 (A): Ensemble spread of the plume area of ENS_REF in km2 (black line, left axis). Ratio between the ensemble spread and the plume area from GOT_REF (red line, right axis). (B): plume area (in km 2 ) from GOT_REF (red line, left axis). Plume area from all the members of the ensemble (black line, left axis). Daily runoff from the Red River and southern rivers (blue line, right axis in m 3 /s).

  Figure V.15 shows that the MAX clusters are equal to the reference clusters (REF on the figure), except on 21 July when the most likely cluster changes from 4 to 3. As a consequence, we estimate that within our experimental conditions the findings of Nguyen-Duy et al. (

Figure V. 15

 15 Figure V.15 Classification of the daily plume patterns in June-July 2015 across the ENS_REF ensemble, using the reference clusters from Nguyen-Duy et al. (2021) and Chapter IV. REF shows the same classification for the GOT_REF run. MAX is the cluster number with the MAXimum likelihood, defined as the most frequent cluster draw within the 50 samples shown underneath.

Figure V. 16

 16 Figure V.16 (A): daily mean surface concentration (unit.m -3 ), computed from 5 -10 July for member 33. (B): same as (A) but for surface wind (m.s -1 ). (C): same as (A) but for surface current (cm.s -1 ). (D, E, F): same as (A, B, C) but for member 34.

Figure V. 17 :

 17 Figure V.17: same as Figure V.16 but computed from 15 -20 July.

  V.15), and the same time of peaks in the surface current spread (Figure V.7A) and in the SSS spread (Figure V.12A). In Section V.3.3 we analyze in more details the response to the wind perturbations in an attempt to explain such a variability.

Figure

  Figure V.19 (A): Relative surface vorticity (10 -6 s -1 ) (contour fill, left color bar) and plume thickness (m) (line, right color bar) on 24 July 2015 of GOT_REF. (B, C, D): same as (A) but for members 18, 38, and 22 of the ENS_REF.

  Figures V.19 B, C, D show the relative surface vorticity and plume thickness on the same day but from 3 members of the ENS_REF. Under different wind conditions, the plume can be either thicker (Figure V.19B), thinner (Figure V.19C), or mostly disappear (Figure V.19D).

From 25 -

 25 30 July, both the wind stress of GOT_REF and the spread of wind stress is high (Figure V.5B). It is also the time when the ratio between the spread of the plume area and the unperturbed plume area is high (Figure V.14A).

Figure V. 20

 20 Figure V.20 (A): The spread of the SSS on 27 July 2015. (B): the quantile of plume occurrence on the same day. The contour lines show the area where the plume occurs in 5 (blue), 25 (black) and 45 (red) times over 50 members of ENS_REF. Cyan color shows the plume area from GOT_REF.

Figure V .

 V Figure V.20A shows the ensemble spread of SSS on 27 July 2015. The spread is highest near the RR mouth (>1) as expected, but it is also high in the northeast of the gulf which can be due to the runoff from northern rivers. The spread is high offshore, but is low near the coast from 19.0 -19.8°N, indicating that the wind may be upwelling favorable in most members.The quantile of plume occurrence is shown in Figure V.20B. Along the coast from 20 -21°N, the plume seems to be stable across the ensemble, as demonstrated by the small difference between the 5-member and 45-member plume occurrence. Along the latitude of 20N, the plume spread is large. For half of the members, the plume can extend to 106.3°E, while in 5 members it can reach around 107.4°E, which is more than 100km from the coast. A large spread is also observed near 19°N.

Figure V. 21

 21 Figure V.21 Surface tracer concentration of GOT_REF (A), the member with the largest plume area (B), the member with the smallest surface plume area (C) on 27 July (arbitrary unit/m 3 ). (D, E, F): Daily mean wind conditions (m/s). (G, H, I): Vertical diffusivity at the second layer from the surface (in m 2 /s). (J, K, L): Vertical profiles of tracer concentration at 20°N (arbitrary unit/m 3 ). (M, N, O): Surface current (cm/s).

Figure

  The wind patterns are similar, which is southeasterly in the south and southwesterly in the north of the gulf. Under upwelling favorable wind, the change of wind speed can have 2 types of impacts.

Figure V. 22

 22 Figure V.22 (A): The spread of the SSS on 20 August 2015. (B): the quantile of the plume occurrence on the same day. The contour lines show the area where the plume occurs in 5 (blue), 25 (black) and 45 (red) times over 50 members of ENS_REF. Cyan color shows the plume area from GOT_REF.

Figure V. 23

 23 Figure V.23 A, B, Cshows the surface tracer concentration of the GOT_REF, the member with the largest plume area (member 40) and the one with the smallest plume area (member 39) on 20 August.In member 39, the small plume area is explained by the strongest wind, around 6 -7 m/s near the coast between 20°N and 21°N (Figure V.23F), similar to the event of the 27 July, and also due to the strong downcoast current (Figure V.23I).In GOT_REF, the coastal wind is the weakest (1 -3 m/s) between 20°N and 21°N (Figure V.23D), however, the plume area does not reach the maximum like for the event on 27 July. Member 40, which has the in-between wind speed (3 -5 m/s, Figure V.23E), has the largest plume area. The largest plume area in member 40 can be explained by the strong offshore transport current between 20 -21°N (Figure V.23H). It means that the runoff plays an important role on the spread of the plume. On July 27, when the runoff is low, the increase of the current due to the wind perturbation can dilute the plume, leading to a smaller plume area. However, on August 20, when the runoff is high, the increase of the current can help the plume to be advected more offshore and increase the plume area.

Figure V. 23

 23 Figure V.23 Surface tracer concentration of GOT_REF (A), the member with largest plume area (B), the member with smallest surface plume area (C) on 27 July (arbitrary unit/m 3 ). (D, E, F): Vertical profiles of tracer concentration at 20N (arbitrary unit/m 3 ). (G, H, I): daily mean wind conditions (m/s). (J, K, L): surface current (cm/s). Please note that the color scales for D, E, F are different.

Figure V. 24

 24 Figure V.24 Example verification rank histograms for ensembles with 8 members. Perfect rank uniformity (also called histogram flatness) is indicated by the dashed lines. Adapted from Wilks (2005).

Figure V. 25 :

 25 Figure V.25: (A): Area to compute the rank histogram for SST (blue box). (B) Area to compute the rank histogram for surface currents (yellow-filled polygon).

First

  , the inconsistency happens when the runoff reaches its peak (Figure V.14). Second, in the model configuration, the maximum temperature of the river runoff is 29°C (see section II.2.2.2), while the SST in summer can reaches over 30°C on average (Figure III.2).

Figure V. 26

 26 Figure V.26 Rank histograms of ENS_REF SST versus OSTIA data over the box, using 51 ranks (Rank 1 and 51 are open rank). F defines the histogram flatness. The y-axis shows the rank. Gray (pink) colorshows the pdf that is lower (higher) than the color scale.

Figure V. 27 :

 27 Figure V.27: Rank histograms of ENS_REF versus HF data, for U component (A) and V component (B). The y-axis shows the rank. There is no data on 20 and 21 June.

Figure V. 27 :

 27 Figure V.27: (A): Spatial mean of ensemble spread (red line) of ENS_REF and spatial mean of RMSE between HF data and ensemble mean of ENS_REF, for U component. (B): same as (A) but for V component. Units are cm/s

3.

  Les régimes du panache du FR ont été identifiés et liés aux processus physiques à l'aide d'une analyse en clusters K-means (KMA).4. Nous avons réalisé un ensemble de simulations pour évaluer la réponse du modèle aux perturbations du forçage du vent. Nous avons analysé la variabilité des incertitudes résultantes dans l'espace et le temps des variables océaniques. L'approche de clustering K-means a été utilisée sur les membres de l'ensemble, donnant accès à une évaluation de la robustesse des résultats du point 3. ci-dessus.

  2.8 𝐾 " = #2𝐸 # 𝑙 # 𝑆 "

	II.1.3 Boundary conditions
	II.1.3.1 Surface boundary
	At the surface, the vertical flux of turbulent kinetic energy (𝐹 % = 𝐾 #	./ &
														Eq. 2.9
	with:
	𝐸 ' : Turbulent Kinetic Energy (TKE)
	𝑙 ' : mixing length
	𝑆 & , 𝑆 # : stability functions (Canuto et al. 2001)
	𝐸 ' and 𝑙 ' are calculated by solving the following equations:
	𝜕𝐸 * 𝜕𝑡	+	𝜕𝑢𝐸 * 𝜕𝑥	+	𝜕𝑣𝐸 * 𝜕𝑦	+	𝜕𝑤𝐸 * 𝜕𝑧	=	𝜕 𝜕𝑧	(𝐾 #	𝜕𝐸 * 𝜕𝑧	) + 𝑃 + 𝐵 -𝜀	Eq. 2.10
			𝜕𝜀 𝜕𝑡	=	𝜕 𝜕𝑧	(	𝐾 # 𝜎 *	𝜕𝐸 * 𝜕𝑧	) +	𝜀 𝐸 *	(𝑐 + 𝑃 + 𝑐 , 𝐵 -𝑐 " 𝜀)	Eq. 2.11
	With 𝜀 is the dissipation rate of the turbulent kinetic energy while P is the production of
	turbulent kinetic energy and B the buoyancy flux. c0 = 0.5544, c1=1.44 , c2 = 1.92. c3 = 1 if
	B ≥ 0 and c3 = -0.52 otherwise (Warner et al., 2005). 𝜎 ' = 1.3 is the Schmidt number for
	dissipation.
	The relationship between the turbulent length, 𝑙 ' , 𝐸 ' and 𝜀, writes:
										𝑙 * =	𝜀 ,/" , 𝐸 * 𝑐 !	Eq. 2.12
	P and B are calculated as:
							𝑃 = 𝐾 # JK	𝜕𝑢 𝜕𝑧	L	"	+ K	𝜕𝑧 𝜕𝑣	L	"	M	Eq. 2.13
									𝐵 =	𝑔 𝜌	𝐾 '	𝜕𝑧 𝜕𝜌	Eq. 2.14
	Eq. 2.13 expresses that the production of TKE is due to the vertical shear of the horizontal
	current.

  (𝜆, 𝜑)𝑐𝑜𝑠 ^𝜔𝑘 (𝑡 -𝑡 0 ) + 𝑉 0,𝑘 + 𝑢 𝑘 -𝐺 𝜂,𝑘 (𝜆, 𝜑)_ Eq. 2.28a𝑓 & 𝑢 ',& (𝜆, 𝜑)𝑐𝑜𝑠 5𝜔 & (𝑡 -𝑡 ' ) + 𝑉 ',& + 𝑢 &

	𝑢(𝜆, 𝜑, 𝑡) = 4
		Eq. 2.28b
	-𝐺 (,& (𝜆, 𝜑)6
	𝑣(𝜆, 𝜑, 𝑡) = 4	𝑓 & 𝑣 ',& (𝜆, 𝜑)𝑐𝑜𝑠 5𝜔 & (𝑡 -𝑡 ' ) + 𝑉 ',& + 𝑢 &
	&	Eq. 2.28c
	-𝐺 ),& (𝜆, 𝜑)6

&

II.2.2 River configuration

	both components of the current are detided based on an online harmonic analysis (as described
	in II.1.4).
	Table II.1: General characteristics of the reference simulation (GOT_REF)
	Name	Features
	Grid points	429 * 600
	Grid size	300m (near the coast), 4500m (near the boundary)
	No. of vertical layers 20
	Boundary conditions	global OGCM
	Atmospheric	ECMWF analysis (1/8°, 3h)
	conditions
	Tidal forcing	FES2014 (K1, O1, P1, Q1, K2, M2, N2, S2, M4)
	Light attenuation	spatially variable (~1m near the coast and 15m offshore)
	Bottom roughness	1 mm
	Timestep	2 minutes
	Period of study	Jan 1, 2010 -Dec 31, 2016
	II.2.2.1 River runoff data
	In the GOT, there are 3 main sources of river runoff: from RR, from other 3 rivers (Ma, Yen,
	Lam) at the south of the Red River delta (hereafter referred to as the 'southern rivers', or SR) , 2.2 is dedicated to its and from 6 rivers at the north of the Red River delta (hereafter referred to as the 'northern description. rivers', or NR). The location of these rivers are shown in Figure II.3.
	The model is run from 2010 to 2016, starting from the ocean state condition as provided by the
	global OGCM on 01/01/2010. The time step is set at 2 minutes. All analyses in this thesis are
	conducted during 2011-2016 (i.e., following a 1-year spinup). The model outputs include daily
	averaged variables as well as instantaneous fields every 12 hours. Unless otherwise specified,

Table II . 2 :

 II2 List of simulations

	Name	Features	Computing specifications
	GOT_REF	Reference configuration	Each simulation is run with 168
	GOT_23	Same as GOT_REF but with L2 = 23m	CPUs on NUWA. CPU time to run GOT_REF: 21 000 hours
	GOT_NORIV	Same as GOT_REF but there is no river.	Disk space for GOT_REF: 2 500 GB.
			Similar specs for GOT_23,
	GOT_NOTIDE	Same as GOT_REF but without tide.	GOT_NORIV, GOT_NOTIDE
			Each simulation is run with 180
	ENS_REF ENS_0.3	2 sets of the ensemble simulations. Details will be presented in Chapter 5.	CPUs on CALMIP. CPU time to run ENS_REF: 75 000h Disk space for ENS_REF: 8 600 GB.

Table III

 III

		Unit in °C
	Data source Configuration	Mean	Min bias	Max bias	Mean R
		absolute bias
	GOT_REF	0.61	-1.28	2.25	0.71
	VITEL data
	GOT_23	0.82	-1.14	3.14	0.57
	GOT_REF	0.81	-1.02	1.34	0.74
	CFO data
	GOT_23	0.89	-0.57	1.78	0.57

.

1 Comparison of GOT_REF and GOT_23 with VITEL and CFO temperature data.

TABLE 1 |

 1 General characteristics of GOT_REF.

	Name	Features
	Grid points	429 * 600
	Grid size	300 m (near the coast), 4500 m (near the boundary)
	No. of vertical layers	20
	Boundary conditions	global OGCM
	Atmospheric conditions	ECMWF analysis (1/8 • , 3 h)
	Tidal forcing	FES2014 (K1, O1, P1, Q1, K2, M2, N2, S2, M4)
	Light attenuation	Varying by space (∼1 m near the coast and 15 m
		offshore)
	Bottom roughness	1 mm
	Timestep	2 min
	Time	2010 -2016

  Wang et al. (2014) for the Copper River in the Gulf of Alaska]

	fmars-08-772139	November 11, 2021	Time: 12:13	# 13
	Nguyen-Duy et al.		Red River Plume Variability

Table V .

 V

	component, the absolute value of unperturbed wind stress is less than 0.02 Nm
		1: List of simulations
	Ensemble simulation	Number of members	𝜀
	GOT_REF	1	0
	ENS_REF	50	0.6
	ENS_0.3	50	0.3

  (Figure V.21 G, H, I). It reaches its highest (lowest) value in the strongest (weakest) wind condition. The vertical profiles of the tracer concentration confirm that the tracer can mix deeper in member 42 (which has strongest wind) to more than 30m, while in member 41 (which has weakest wind) the tracer is mixed in less than 15m (Figure V.21 J, L).

  De septembre à mars, le panache est étroit, étendu le long de la côte et 'attaché au fond' (bottomattached) en raison des vents de downwelling de la mousson d'hiver. A partir de mai, la mousson se transforme en vent de sud-ouest, poussant le panache vers le nord et le détachant de la côte et du fond. Dans ce cas, le panache est advecté en surface avec une épaisseur de 2-3 à 7-8 m. Nous observons que les événements où le panache se détache de la côte à 19°N coïncident avec la température à la côte plus froide qu'en mer. Ceci soutient l'hypothèse que les vents de mousson du sud génèrent localement un upwelling côtier. Parfois, en mai-août, le panache s'épaissit au large jusqu'à 12-15 m. Cet approfondissement peut être dû à des vents forts qui augmentent le mélange au large dans la région d'upwelling ou à l'augmentation du débit des rivières.A la mi-août, nous observons un fort approfondissement du panache sur le plateau (106-106.5°E). Ces événements d'approfondissement semblent être liés à un tourbillon anticyclonique saisonnier récurrent se développant près de 19°N. Dans notre simulation, ce tourbillon se produit généralement en août, lorsque la direction du vent est du sud-ouest, ce qui est favorable à l'upwelling, et que le débit est élevé. Ensuite, ce tourbillon disparaît lorsque le courant côtier vers le sud se développe à nouveau et/ou que le vent n'est plus upwellingfavorable. Le mécanisme de formation de ce tourbillon n'est toujours pas clair. Dans GOT_NORIV, des gradients de vorticité peuvent être représentés, suggérant une faible structure anticyclonique, d'amplitude beaucoup plus faible que dans la simulation de référence et avec une durée de vie plus courte. L'anticyclone semble fortement connecté au panache de la rivière Lam. La topographie du fond peut être un autre facteur de forçage puisque le tourbillon est situé à proximité d'une île.

https://sextant.ifremer.fr/Donnees/Catalogue#/metadata/2b192671-1129-409c-adae-e2335759f290

Frontiers in Marine Science | www.frontiersin.org

November 2021 | Volume 8 | Article 772139

November 2021 | Volume 8 | Article

https://numpy.org

ACKNOWLEDGMENTS

CFO CTD data used in this study was collected by the Vietnam Center for Oceanography (CFO, VASI) in the framework of the United States-Vietnam collaboration project: "Gulf of Tonkin Circulation study" (NICOP N62909-15-1-2018), financed by the Office of Naval Research (ONR). The VITEL (VIetnam TELédétection remote sensing) campaign has been funded by CNES; we are grateful to S. Ouillon for providing the T/S data. The historical tidal gauges dataset stems from the International Hydrographic Organization (https://www.iho.int/, last access: November 2, 2019) and is available upon request at https://www.admiralty.co.uk/ukho/tidal-harmonics (last access: November 2, 2019). The altimetric tidal constituents were provided to us by F. Lyard (LEGOS). We also acknowledge F. Lyard for the use of his code to compute the Rossby Radius and for his expertise on the model's bathymetry field. This study has been conducted using ECMWF product (www.ecmwf.int) and E.U. Copernicus Marine Service Information (product: GLOBAL_ANALYSIS_FORECAST_PHY_001_024, available at https://resources.marine.copernicus.eu/?option=com_csw& view=details&product_id=GLOBAL_ANALYSIS_FORECAST_ PHY_001_024). Many thanks to C. Nguyen and D. Gazen for the help with the model and computing facilities at LAERO (Toulouse). The SYMPHONIE code and the code to compute the Rossby radius are distributed by the SIROCCO French national service website (https://sirocco.obs-mip.fr/), hosted by the Observatoire Midi-Pyrénées, Toulouse, France.

FUNDING

TuN-D Ph.D. thesis is funded by the "Programme d'Excellence" of the French Embassy of Vietnam and by IRD ARTS program. This work is a contribution to the LMI LOTUS project (http: //lotus.usth.edu.vn) and to the SWOT COCTO-FO project funded by CNES (French Spatial Agency). This work has used HPC resources from CALMIP/France (grant 2020-p1119) and HILO at the University of Science and Technology of Hanoi (USTH). Support from IRD, CNES, CNRS, and USTH is gratefully acknowledged.

DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article (i.e., the model simulations) will be made available by the authors without undue reservation.

AUTHOR CONTRIBUTIONS

TuN-D, NA, PM, and PDM-F contributed to the conception and scientific design of the study. PM provided the model. VP, MH, TD, and PM developed the previous GOT configuration used for this analysis. MT provided the HFR data. TuN-D, NA, and FT performed most of the analyses. TuN-D wrote the first draft of the manuscript. NA wrote sections of the manuscript. All authors actively contributed to manuscript revision, read, and approved the submitted version.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2021.772139/full#supplementary-material

Conflict of Interest:

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.

Publisher's Note: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

APPENDICES APPENDIX A: POST-PROCESSING OF THE CMEMS T,S FIELDS

On the continental shelves, the turbulence generated by the friction of the tidal currents on the bottom is generally significant. Near the bottom, the temperature and salinity stratification are in principle homogenized over a layer that can be up to several tens of meters, highly depending on the strength of the tides. CMEMS fields provide boundary conditions for temperature, salinity, surface elevation and currents associated with non-tidal processes. Because tides are not taken into account in the CMEMS simulation, the stratification is unfortunately insufficiently mixed near the bottom. As the errors induced on T and S are potentially important, we make a modification to the CMEMS T,S fields before introducing them at the open boundaries of our model. This modification consists in homogenizing T and S near the bottom, over a thickness deduced from the expected theoretical height of the tidally induced turbulent bottom layer. This height is estimated from the friction velocity associated with the tidal current, i.e. 𝐷 = 0.1𝑢 * /𝑓 (Csanady, 1982) where 𝑓 is the Coriolis parameter. The friction velocity, 𝑢 * , is calculated from the tidal harmonics of the depth-averaged current provided by FES2014, and a drag coefficient of 1e -3 .

APPENDIX B: RIVER BOUNDARY CONDITIONS

In this appendix, we describe the way the inflow of river water is implemented at the most upstream point of the channel (if any) or directly at the river mouth into the ocean in case there is no channel. To simplify, we will call it the "river boundary condition at the river input point". (Note that for a given river, in case there is no channel, there can be several input points to represent a river mouth wider than the actual grid size; in that case, the boundary condition as described below is easily adapted.)

At the river input point, the river boundary current is consistent with the runoff of the river as in Reffray et al (2004), while the vertical shape of the current is given by a gradient condition allowing the current to adjust by itself to the dynamics of the river downstream (rather than arbitrarily prescribing a possibly inappropriate vertical profile). The boundary condition is applied in two steps. Step 1 is given by the above-mentioned gradient condition, leading to a provisional solution of the current:

if current in downstream direction of the river 𝑢 * = 0 otherwise where x is the direction of the river flow and u* a provisional solution of the current.

Step 2 ensures that the total flow is consistent with the expected runoff from the river. In practice, the river boundary current is proportional to 𝑢 * , i.e. 𝑢 = 𝛼 𝑢 * , where 𝛼 is a constant used to satisfy the constraint on the prescribed river runoff:

where 𝛥𝑦 is the cellbox width at the river input point, h the bottom depth, 𝜂 the sea surface height and F the river runoff.

Appendix C: Equation of the distance and procedure to identify the number of clusters.

The distance from a member (X) to the centroid (A) is calculated using Euclidean distance. The distance d(X,A) is computed for a given variable v defined at every model grid point i (i=1,N).The equation of the Euclidean distance is:

where Xi: value of variable v at grid point i of the member X Ai: value of variable v at grid point i of the centroid of cluster A N: the total number of grid points

As explained in section 2.4, the variable v at point i takes only two values: Xi=1 if the plume is present, Xi=0 otherwise.

In the K-means method, the number of clusters has to be chosen a priori. The number is usually a trade-off between the need to reduce the dimension of the ensemble of scenes to be analyzed and the search for meaningful distinct clusters. The number of clusters can be set subjectively (e.g. Solabarieta et al., 2015), based on a priori knowledge of the studied process (here the plume variability), or can rely on information criteria (Sonnewald et al., 2019).

In this study, since KMA is based on the distance to the cluster centroid, in order to evaluate how many clusters are suitable, we compute the sum of all the distances between the members and the centroids. The more clusters there are, the smaller the total distance will be decreased.

We applied the KMA for 8 cases, from 2 to 9 clusters. The total distance and the distance reduction are shown in Figure C1. We can clearly see that the passage from 2 to 3 clusters induces a significant reduction of the distance (more than 2x10 6 ). Similarly, from 3 to 4 clusters, the distance reduction is still large. However, from 4 to 5 clusters, the distance reduction is decreased significantly (~ 0.8x10 6 ). It means that using 4 clusters is suitable in terms of the total distance. simulations, to compute statistics from the ensembles, such as spread, quantiles, and to analyze the probabilistic distribution. Ensemble methods in oceanography were introduced about three decades ago for data assimilation purposes (e.g. Evensen 1994). Since then, they have been used as well for uncertainties estimates or budgets in open ocean (e.g. [START_REF] Andreu-Burillo | Assessing oceanmodel sensitivity to wind forcing uncertainties[END_REF][START_REF] Lucas | Stochastic study of the temperature response of the upper ocean to uncertainties in the atmospheric forcing in an Atlantic OGCM[END_REF] or coastal ocean studies (e.g. [START_REF] Auclair | Space-time structure and dynamics of the forecast error in a coastal circulation model of the Gulf of Lions[END_REF][START_REF] Kim | Toward an uncertainty budget for a coastal ocean model[END_REF], and generalized to coupled models such as physical-biogeochemical ones (e.g. Vervatis et al. 2021a). Probabilistic modeling is also implemented in numerical weather prediction systems, such as ECMWF, while ensembles of models are used to assess seasonal to interannual predictions (e.g. [START_REF] Kirtman | The North American Multimodel Ensemble: Phase-1 Seasonal-to-Interannual Prediction; Phase-2 toward Developing Intraseasonal Prediction[END_REF], Herrmann et al. 2021).

Supplementary materials

As said before, the work of [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF] in Chapter IV has revealed the role of wind on the variability of the plume. However, until now, the robustness of these results has not yet been assessed. In particular, since the plume variability is strongly influenced by the wind, as evidenced by the cluster analysis, the impact of uncertainties in the surface forcing may be critical. More generally, the coastal and shelf sea dynamics respond mainly to the wind action, sometimes combined with the effect of tides and river runoff. For instance, by perturbing several sources of errors (atmospheric forcing, model parameterizations), [START_REF] Vervatis | Assessment of a regional physical-biogeochemical stochastic ocean model. Part 1: Ensemble generation[END_REF] suggest that in the Bay of Biscay, the main source of error is the wind forcing. This is also the conclusion of other studies, for example [START_REF] Jordà | Characterization of error dynamics in a 3D coastal model of the Catalan sea using stochastic modelling[END_REF] in the Catalan sea or Auclair et al. (2002) for the Gulf of Lions. Therefore, this study aims at assessing the robustness of the Chapter IV results and conclusions with respect to an important source of uncertainties, namely the wind forcing. This in turn translates into the following three questions:

1. What is the space and time variability of the ensemble spread of the surface ocean variables ? In particular, does the model sensitivity to wind uncertainties depend on the runoff regime ?

2. Can we identify the physical processes at work to explain the sensitivity and its variability ?

3. Is the cluster analysis of [START_REF] Nguyen-Duy | Variability of the Red River Plume in the Gulf of Tonkin as Revealed by Numerical Modeling and Clustering Analysis[END_REF] robust with respect to the simulated wind uncertainties ?

Besides, the ensemble approach raises several methodological questions, three of which we address in this study:

• How to generate wind perturbations that are indeed representative of ECMWF wind errors ? • How (by means of which metrics) do we define the sensitivity to wind uncertainties with an ensemble? • Is the generated ocean ensemble significantly representative of model errors ?

In order to do that, we will run an ensemble of simulations and assess the model response to perturbations added to the wind forcing. This category of methods has been used by several authors [START_REF] Ayoub | Estimating uncertainties on a Gulf Stream mixed-layer heat budget from stochastic modeling[END_REF][START_REF] Ghantous | Ensemble downscaling of a regional ocean model[END_REF][START_REF] Jordà | Characterization of error dynamics in a 3D coastal model of the Catalan sea using stochastic modelling[END_REF], Quattrocchi et al. 2014[START_REF] Vandenbulcke | A stochastic operational forecasting system of the Black Sea: Technique and validation[END_REF], among many others). Firstly, we will try to assess the uncertainty of the wind forcing. In some studies (Quattrocchi et al. 2014, Vervatis et al. 2016[START_REF] Ghantous | Ensemble downscaling of a regional ocean model[END_REF], the perturbations of the wind forcing are based on the author's assumption on the

V.3 Results

In this section, we will assess the impact of the wind perturbation to the ocean. The assessments are mostly based on the ensemble spread, which is computed as the standard deviation along the ensemble with respect to the ensemble mean. 2021) (the base of the plume is the depth where the concentration is equal to 7 units/m 3 ). In June, the plume is about 2 m deep and detached from the coast. Then, the thickness increases and the plume extends further offshore, reaching the maximum thickness of around 10m on 20-25 July. In August, the plume thickness varies from 2 to 7m.

V.3.1 Ensemble spread of the standard model variables

The spread of the plume thickness is shown in Figure V.18B. Due to the wind perturbation, the plume is not only changed vertically but also changed horizontally. Depending on the perturbation, the plume can stick nearer to the coast or extend further offshore. Due to the horizontal displacement of the plume, it is complicated to assess the spread of the ensemble thickness. After August 20, east of 106.3°N, the spread is nearly equal to 0 not because the plume is stable across the ensemble, but due to the fact that only the plume of 1-2 members can reach that longitude (Figure V.18C). It means that the presence of the plume in these points is actually unstable across the ensemble, which is consistent with the high spread of the plume area analysis above.

VI.1 Conclusion

Gulf of Tonkin is a shallow region (depth <60m) of the southwest Pacific, under the influence of the Red River (hereafter RR) discharge, of a strong atmospheric variability (monsoons, typhoons) and of tides. The Red River delta has a total population of 22.5 millions and an average population density of 1060 inhabitants/km 2 , which is nearly four times higher than the national average for Vietnam. The very high population density exerts a strong pressure on the rivers and their environment.

This work is part of the LOTUS project. LOTUS is a joint Vietnam-France project, which aims to better understand the coastal ocean dynamics under the influence of atmospheric variability and of the delta. In the PhD thesis of V. [START_REF] Piton | Du Fleuve Rouge au Golfe du Tonkin: dynamique et transport sédimentaire le long du continuum estuaire-zone côtière[END_REF] at LEGOS (Toulouse, France), a highresolution configuration in the Gulf of Tonkin and downstream part of the Red River delta has been set-up to study tides with the 2D spectral TUGO model and basin-scale circulation with the 3D SYMPHONIE model. In this study, the main objective was to provide a better understanding of the variability of the RR plume in the Gulf of Tonkin at different scales using an updated version of the configuration with SYMPHONIE.

The main scientific questions which I attempted to address are the following:

1. How can high-resolution numerical simulations reproduce the Gulf of Tonkin dynamics and the Red River plume variability at different scales in the study years? 2. Which criteria can be used to identify the RR plume in the GoT from the model outputs? 3. How can we formally characterize the sequence of regimes which the RR plume goes through in the study years? Can we link those regimes to specific processes and forcings ? 4. How can we assess the sensitivity of model variables and in particular on the variability of the river plume, to the main sources of uncertainty, in particular the atmospheric forcing uncertainties?

To address those questions, we adopted the following tool-based approaches:

1. We set up and run an upgraded version of the SYMPHONIE model configuration created as part of the thesis of V. Piton. A reference simulation (GOT_REF) is run with realistic forcing conditions over a 7-year period (2010)(2011)(2012)(2013)(2014)(2015)(2016). Additionally, a twin simulation without river forcing (GOT_NORIV) is performed to assess the impact of the river runoff on the coastal circulation and a twin simulation without tides (GOT_NOTIDE) is run to evaluate the impact of tides on the main patterns of the river plume variability. 2. The river plume was identified using different methods, allowing us to select the most appropriate one given our purposes. The selected method is based on passive tracers injected in the rivers. 3. The RR plume regimes were identified and linked to physical processes with the help of K-means clustering analysis (KMA). Pour répondre à ces questions, nous avons adopté les approches suivantes basées sur différents outils :

CONCLUSION GÉNÉRALE

1. Nous avons mis mis à jour et amélioré la configuration de SYMPHONIE créée dans le cadre de la thèse de V. Piton. Une simulation de référence (GOT_REF) est produite dans des conditions de forçage réalistes sur une période de 7 ans (2010)(2011)(2012)(2013)(2014)(2015)(2016). En outre, une simulation jumelle sans forçage fluvial (GOT_NORIV) est mise en place pour évaluer l'impact de l'apport d'eau des rivières sur la circulation côtière et une simulation jumelle