Manh Cuong

Alexei Sentchev

Stefano Berti

Nadia K Ayoub

Tung Nguyen-Duy

Kim Nguyen

Cuong

Assessment of relative dispersion in the Gulf of Tonkin using numerical modeling and HF radar observations of surface currents

Keywords: Relative dispersion, Scale of motion, HF radar measurements, Numerical modeling, Gulf of Tokin

Particle pair statistics from synthetic drifter trajectories reconstructed from realistic, high-resolution numerical simulations (SYMPHONIE model) and HF radar velocity measurements are used to investigate the dispersion properties in the Gulf of Tonkin (GoT). This study takes an approach based on two-particles statistics providing the relative dispersion, relative diffusivity and Finite Size Lyapunov Exponent (FSLE) estimates. In the GoT, the relative dispersion follows the predictions from the theory of two-dimensional turbulence with two inertial subranges identified in the kinetic energy spectrum with the spectral slopes -5/3 and -3. The time evolution of dispersion shows an exponential growth during 5 8 days, followed by a power law regime during the next 6 to 20 days. Fixed-length indicators from the relative diffusivity and the FSLE reveal a local dispersion at large and intermediate scales (above Rossby radius of deformation) and non-local dispersion in sub-mesoscale range (below Rossby radius of deformation). The effect of river runoff on the local hydrodynamics and dispersion processes is assessed using the numerical model simulations without river discharge. The results show that in the model, the river plume, when present, highly impacts the Lagrangian statistics. High gradients of buoyancy reinforce the sub-mesoscale circulation in a large region along the Vietnamese coast and modify the scales and intensity of turbulent dispersion. However, a clear change of dispersion regime in the sub-mesoscale range is not identified, suggesting that the mesoscale circulation in the GoT largely governs particle spreading even at small scale.

Introduction

Transport and dispersion processes in the ocean have practical importance since they play a major role in the functioning of marine ecosystems by carrying physical and chemical quantities such as heat, salt, nutrients, as well as biological relevant tracers, matter and marine debris. For example, in the North Sea, the transport of the flatfish larvae from the spawning areas to the nurseries and mortality rate are heavily impacted by the hydrodynamics processes and the environmental factors [START_REF] Lacroix | Modelling dispersal dynamics of the early life stages of a marine flatfish (Solea solea L.)[END_REF]. In the coastal zone, the knowledge of transport properties is even more important for assessing the dispersal of anthropogenic pollutant tracers released accidentally [START_REF] Corrado | General characteristics of relative dispersion in the ocean[END_REF]. Mapping the transport pathways of materials advected by ocean currents is of great interest for understanding the interconnection between geographical regions [START_REF] Haza | Transport properties in small-scale coastal flows: relative dispersion from VHF radar measurements in the Gulf of La Spezia[END_REF]van Sebille et al. 2011van Sebille et al. , 2015)), ph sical and relevant biogeochemical processes (d Ovidio et al., 2015). However, prediction of the Lagrangian transport remains challenging since the oceanic flow is turbulent in a wide range of spatial and temporal scales. The oceanic coherent structures, which we refer to as jets, eddies, fronts, arise from multi-scale interactions of physical processes in the ocean and at its interfaces with the atmosphere and the continent. The instabilities of these structures, due to the flow turbulence or episodic forcing (e.g., wind), result in chaotic motions and may dictate the transport statistics [START_REF] Haza | Impact of submesoscales on surface material distribution in a gulf of Mexico mesoscale eddy[END_REF].

Lagrangian framework is naturally linked to the investigation of transport since the Lagrangian observations track the fluid parcels and sample the spatiotemporal evolution of their properties rather than measuring them at a fixed point in the Eulerian approach [START_REF] Essink | Can we detect submesoscale motions in drifter pair dispersion[END_REF]. Trajectory maps can provide the details of the circulation and dispersion of materials in the flow field. Since an individual trajectory in the turbulent flow field is largely unpredictable and depends strongly on the initial conditions, the ensemble average over many trajectories is needed to provide the statistical description of the flow field (La Casce et al., 2008;[START_REF] Essink | Can we detect submesoscale motions in drifter pair dispersion[END_REF].

Characteristics of the oceanic turbulent flow based on the Lagrangian approach have been described in many studies, e.g., in the North Atlantic (Lumpkin and Elipot, 2010), in the South Atlantic, [START_REF] Berti | Lagrangian Drifter Dispersion in the Southwestern Atlantic Ocean[END_REF], in the Nordic Seas (Kos alka et al., 2009), in the Gulf of Me ico (LaCasce and Ohlmann, 2003;[START_REF] Asaro | Ocean convergence and the dispersion of flotsam[END_REF], and in the Mediterranean Sea [START_REF] Haza | Transport properties in small-scale coastal flows: relative dispersion from VHF radar measurements in the Gulf of La Spezia[END_REF][START_REF] Schroeder | Relative dispersion in the Liguro-Provençal basin: From sub-mesoscale to mesoscale[END_REF][START_REF] Schroeder | Targeted Lagrangian sampling of submesoscale dispersion at a coastal frontal zone[END_REF]. These studies aimed to explore how pollutants or biological tracers disperse in the marine environment and how the energy is transferred across scales of motion. The relative dispersion, which might be estimated from the mean squared separation of particle pairs, is used to characterize the underlying turbulent flow field. Among the Lagrangian statistics, pair dispersion gains importance since the spreading rate of pairs of particles is sensitive to the amount of energy distribution within different scales. The spreading of particle pairs depends, in principle, on the velocity difference at the length scale corresponding to the pair separation distance (Lumpkin and Elipot, 2010). Therefore, the statistics based on relative dispersion are useful for the assessment of mixing in the turbulent flow field at small and intermediate scales and can characterize the most important physical processes controlling the underlying fluid motions (Meyerjürgens et al., 2019).

The relationship between the oceanic turbulent flow and the relative dispersion statistics has been extensively discussed (see, e.g., LaCasce, (2008); [START_REF] Foussard | Relative dispersion in generalized two-dimensional turbulence[END_REF]; [START_REF] Berti | Lagrangian pair dispersion in upper-ocean turbulence in the presence of mixed-layer instabilities[END_REF]). It was shown that the dispersion properties can be theoretically predicted from the kinetic energy spectrum, E(k) ~ k -, of the underl ing flow field. Specificall , the relative dispersion can be either local if 1 < < 3, or non-local, in the case of 3. In a number of studies using the relative dispersion statistics [START_REF] Koszalka | Relative dispersion in the Nordic Seas[END_REF]Lumpkin and Elipot 2010;[START_REF] Poje | Submesoscale dispersion in the vicinity of the Deepwater Horizon spill[END_REF]van Sebille et al. 2015;[START_REF] Sansón | Surface relative dispersion in the southwestern Gulf of Mexico[END_REF][START_REF] Beron-Vera | Statistics of simulated and observed pair separations in the Gulf of Mexico[END_REF][START_REF] Corrado | General characteristics of relative dispersion in the ocean[END_REF][START_REF] Essink | Can we detect submesoscale motions in drifter pair dispersion[END_REF], the geographic dependency of the stirring regime is found across different ocean regions. Contradictory results in some regions have been reported revealing a non-local dispersion regime at scales below the Rossby radius of deformation and local dispersion in the same scale range. For example, Meyerjürgens et al. (2019) identified significantly different dispersion regimes in the vicinity of the tidal front in the North Sea. Scaling of the relative dispersion indicated a non-local dispersion regime while the Finite Size Lyapunov Exponent (FSLE) showed a local dispersion. At large spatial scales, the existence of different dispersion regimes was revealed. In some studies (e.g., [START_REF] Schroeder | Relative dispersion in the Liguro-Provençal basin: From sub-mesoscale to mesoscale[END_REF] in the Liguro-Provençal basin), a transition to the diffusive regime was identified, whereas other studies indicated a Richardson-like dispersion, most likely due to the effect of large-scale velocity shear (e.g., LaCasce and Ohlmann 2003 in the Gulf of Mexico). In addition, characterization of the dispersion regime based on the observations from drifting buoys alone is practically difficult and requires many pairs of Lagrangian drifters whose number is always limited in ocean regions. The lack of observations leads to differences in assessment of the dispersion regime in the world oceans (i.e., Haza et al., 2008;[START_REF] Schroeder | Targeted Lagrangian sampling of submesoscale dispersion at a coastal frontal zone[END_REF]. Recently, advances in observations of ocean dynamics by remote sensing techniques such as altimetry, coastal radars, etc., and the improvement in the ocean modeling accuracy have opened a new opportunity for assessment of the oceanic dispersion properties (using synthetic drifters).

The present study is performed in the Gulf of Tonkin (GoT), a shallow, semi-isolated gulf located in the northwestern part of the South China Sea (Fig. 1). As demonstrated in recent studies, the hydrodynamics in the GoT is controlled by the Asian monsoon system, complex bathymetry and exchanges with the South China Sea at the south and through the Hainan Strait (e.g., [START_REF] Tran | Multi-scale variability of circulation in the Gulf of Tonkin from remote sensing of surface currents by high-frequency radars[END_REF][START_REF] Piton | Influence of winds, geostrophy and typhoons on the seasonal variability of the circulation in the Gulf of Tonkin: A high-resolution 3D regional modeling study[END_REF][START_REF] Rogowski | Air-Sea-Land Forcing in the Gulf of Tonkin: Assessing Seasonal Variability Using Modern Tools[END_REF].

The freshwater input from the Red River impacts the hydrology and circulation along the Vietnamese coast with a strong seasonality as shown recently for instance by Nguyen-Duy et al. (2021). In coastal and shallowwater seas, a combination of tidal and wind forcing, in conjunction with complex topography, can lead to the formation of sub-mesoscale eddies (e.g., [START_REF] Zimmerman | Dynamics, diffusion and geomorphological significance of tidal residual eddies[END_REF]. Additionally, large freshwater discharge enhances the density gradients and can also generate sub-mesoscale motions that play an important role in organizing the flow structures and materials transport.

In this work, we employ the realistic, high-resolution numerical model SYMPHONIE and velocity measurements by high-frequency radars to investigate the dispersion properties of the flow field. We adopt an approach based on particle pairs statistics to quantify different scales of the flow field variability and associated dispersion regimes, in the GoT. This article is organized as follows. In section 2, we briefly describe the data and methods used for dispersion quantification. The description of the flow fields during two selected periods is given in section 3. In section 4, results of the assessment of the turbulent dispersion based on the coastal radar velocity measurements and the numerical model simulations are presented. We also assess the effect of river discharge on the Lagrangian statistics. We discuss the results in section 5 and summarize them in section 6.

Materials and methods

HF radar velocity measurements

A coastal radar network, which consists of three radar sites within the GoT since 2012 (Fig. 1). The radar network was configured to operate at a frequency of 4.625 MHz (which was changed to 5.25 MHz in 2018). The radial velocity component is measured by the radars at the effective depth of ~2 m with an instrumental error of ~0.07 m s -1 . The following configuration has been used at each radar site to retrieve the radial velocities of the surface currents: 5.8 km along the beam and azimuthal spacing 5°. The temporal resolution was set to 1 h. The surface HFR total velocities were mapped using EOF/Variational interpolation methods (Yaremchuck and Sentchev et al., 2011). The resulting product consists of hourly maps of the two current components (zonal and meridional) on a regular 6-km grid as described in [START_REF] Tran | Multi-scale variability of circulation in the Gulf of Tonkin from remote sensing of surface currents by high-frequency radars[END_REF]. The data from the HNDU site were not used because of the hardware failures during the study period in 2015.

Numerical model configuration

We use the numerical simulations of Nguyen-Duy et al. ( 2021) made with SYMPHONIE model.

SYMPHONIE is a hydrostatic, numerical model which solves the primitive equations on a curvilinear bipolar Arakawa C-grid with regular sigma vertical levels [START_REF] Marsaleix | Considerations on open boundary conditions for regional and coastal ocean models[END_REF][START_REF] Marsaleix | Energy conservation issues in sigma-coordinate free-surface ocean models[END_REF] surveys in order to characterize the complex topography, with many islands and islets located within the region and especially the coastal area adjacent to the Red River delta (Piton et al., 2019). Intertidal areas are simulated with a wetting and drying scheme. The k-turbulence closure scheme is implemented as in Michaud et al.

( 2012).

Tidal motions at the open boundaries are prescribed from the sea surface elevation and barotropic currents of 9 tidal constituents from the FES2014 atlas [START_REF] Lyard | FES2014 global ocean tide atlas: design and performance[END_REF]. The tidal astronomical potential is also considered at every grid point of the model. In addition, the 1/12° spatial resolution, daily averaged sea surface height (SSH), velocity components temperature, and salinity provided by the Copernicus Marine Systems (CMEMS) global model are also used at the boundaries.

At the sea surface, the boundary conditions include the 3-hourly data for wind, precipitation, solar energy, atmospheric temperature and pressure from the Operational ECMWF product. The fluxes of momentum, heat, and freshwater are computed internally in the model using the bulk formulae of Large and Yeager (2004).

In the reference simulation (hereafter GoT-Ref) the river discharge from 16 rivers is introduced in the model and configured as realistically as possible according to the method described in Nguyen-Duy et al. (2021).

Daily freshwater runoff is used from the main seven distributaries of the Red River. For other rivers, climatological data are used. An additional configuration without the river discharge (GoT-noriv configuration, see Nguyen-Duy et al., 2021) is used to assess the effect of rivers discharge on the dispersion properties.

The model was extensively validated over the period 2011-2016 and the results of the validation are given in [START_REF] Piton | Influence of winds, geostrophy and typhoons on the seasonal variability of the circulation in the Gulf of Tonkin: A high-resolution 3D regional modeling study[END_REF] and Nguyen-Duy et al. (2021). In our study, we performed an additional comparison of the model simulation with the current velocity derived from HFR measurements for two periods (August and November 2015, Fig. 2). The correlation between the current velocity derived from both data sources remains larger than 0.7 with a slightly higher value obtained in November (0.85). A relative difference of about 6-10% quantifies the disagreement between the HFR data and the numerical simulations, for both the latitudinal (u)

and longitudinal (v) velocity components. The model velocities were found to be slightly smaller than those derived from HFR measurements.

Lagrangian diagnostics

In this study, we used the velocity data from the model simulation covering one-year long period in 2015.

Virtual particles organized in clusters were advected using oceanParcels Lagrangian framework (https://oceanparcels.org) by the velocity field from the model at one hour resolution. These Lagrangian tracers, with a total number of 24065 particles, were uniformly distributed in the surface layer at 2 m depth.

Each cluster included one center particle and two satellite pairs along x-(longitude) and y-axis (latitude) with the initial separation 0 = 0.5 km. The clusters were seeded every 2 km across the domain in order to keep the initial pair separation as low as possible and to not exceed a reasonable computational cost. The particles are advected horizontally using the 4 th order Runge-Kutta scheme.

In this study, we focus on assessing the statistics of Lagrangian particle pairs advected in the flow fields derived from the HFR measurements and the model simulation during two seasons (summer and late fall). The Lagrangian quantities analyzed include relative dispersion, relative diffusivity and Finite Size Lyapunov Exponent (FSLE).

Relative dispersion is a common metric for assessing the dispersion of tracers. It is defined as the average, over all particle pairs, of the squared separation between particles i and j in a pair at time t:

𝛿 2 𝑡 1 -1 ∑ ‖𝑥 𝑡 𝑥 𝑡 ‖ 2 , ( 1
)
where N is the number of particles. The pairs were identified using the particles with initial separations below 2 km. The relative dispersion behavior depends on the level of turbulence in the flow field. The latter is reflected in the kinetic energy spectrum which shows a power law distribution as a function of the (horizontal) wavenumber k=|k |: E(k) ≃ k -. At the early stage when the pair separation is small, the particle pair velocities are close and the growth of relative dispersion follows a ballistic behavior, δ 2 ~ t 2 . At larger times, relative dispersion grows exponentially in the case of non-local dispersion, when > 3. Otherwise, if 1 < < 3, relative dispersion is local and scales as a power law δ 2 ≃ 𝑡 3- [START_REF] Foussard | Relative dispersion in generalized two-dimensional turbulence[END_REF]. Richardson super diffusive regime corresponds to δ 2 ~ t 3 where the pair separation is controlled by eddies with a size similar to the separation distance of particles. When the particle separation reaches scales much larger than the largest characteristic flow scale, the pair velocities are considered uncorrelated. Relative dispersion in this case is expected to scale as δ 2 ≃ t, corresponding to a diffusive behavior.

In addition, we diagnose the relative diffusivity κ, which represents the growth rate of the average squared separation. This quantity is defined as the time derivative of relative dispersion,

𝜅 1 2 ⟨𝛿 2 𝑡 ⟩ 𝑡 , (2)
and it can be useful to assess the mixing properties and the corresponding spreading regime. For this purpose, it is often interesting to consider κ as a function of the separation distance δ, rather than of time. Depending on the shape of the given kinetic energy spectrum, κ(δ) scales as δ 2 for the non-local dispersion regime and as δ 4/3 (when =5/3) for the local dispersion one.

The third quantity, the Finite Size Lyapunov Exponent (FSLE), is estimated by measuring the time τ, averaged over all particle pairs, needed to separate particles in a pair from a distance δk to a distance δk+1 = r δk, with r > 1 (k=1, Nk and Nk is the number of pairs chosen in a way to span the relevant range of scales):

𝜆 𝑙𝑛 (3)
Values of FSLE can be used to distinguish a local regime of dispersion from a non-local one [START_REF] Schroeder | Relative dispersion in the Liguro-Provençal basin: From sub-mesoscale to mesoscale[END_REF]. Compared to relative dispersion, the FSLE, being a fixed-length (instead of fixed-time) indicator, better allows to disentangle contributions from turbulent features of different size [START_REF] Artale | Dispersion of passive tracers in closed basins: Beyond the diffusion coefficient[END_REF]. We implemented a method called fastest-crossing for the computation of the FSLE (Lumpkin & Elliot, 2010). A sensitivity study of the results to the choice of r was performed with the fixed number of particles in the area.

It was found that varying r from 1.2 to 1.6 does not change the FSLE slope considerably [START_REF] Poje | Submesoscale dispersion in the vicinity of the Deepwater Horizon spill[END_REF] and the value of r = 1.4 was chosen. The Lyapunov exponent quantifies the exponential rate of separation and is often used for the identification of the nonlocal regime of pair separation [START_REF] Haller | Lagrangian coherent structures and mixing in two-dimensional turbulence[END_REF]. At very small separations, the FSLE is expected to recover such, scale-independent, value. The behavior of the FSLE in a turbulent flow field with spectrum E(k) = k -5/3 typically is as follows: ≃ constant (exponential growth and nonlocal dispersion regime) at very small scales, 𝜆 ≃ 𝛿 -2/3 for local dispersion (Richardson regime) in the inertial range and 𝜆 ≃ 𝛿 -2 for the diffusive regime expected at separations much larger than the largest flow scales [START_REF] Berti | Lagrangian Drifter Dispersion in the Southwestern Atlantic Ocean[END_REF][START_REF] Corrado | General characteristics of relative dispersion in the ocean[END_REF]. A ballistic regime -1 can also be expected in the presence of intense mean currents.

Horizontal stirring from FSLE maps

In addition, the spatial distribution of FSLEs has been obtained for two selected periods representing two seasons (in August and November). The backward FSLE calculations were performed for 30 days with a timestep of 1 hour, following the method proposed by d 'Ovidio et al. (2004). The advantage of this technique is to restrict the analysis of dispersion within relevant structures among selected length scales. The optimum choice for the spatial FSLE amplification factor 𝛿 𝑓 /𝛿 was discussed several times in previous studies [START_REF] Haza | Transport properties in small-scale coastal flows: relative dispersion from VHF radar measurements in the Gulf of La Spezia[END_REF][START_REF] Berta | Surface transport in the Northeastern Adriatic Sea from FSLE analysis of HF radar measurements[END_REF] We compute the exponential growth rate, 𝜆 𝑥 :

𝜆 𝑥 The current vectors are interpolated on a regular grid and plotted every 15 km.

Description of the Eulerian flow fields

Figure 2 shows the 30-days averaged surface current fields from the HF radar measurements and from model simulations in August and November (Fig. 2). During summer, the prevailing wind is from the southwest direction and remains persistent until September 12 when a typhoon hit the region. Interestingly, during August, the surface currents from the model show a pronounced northeast flow in the northern part of the gulf above 19°N which roughly follows the 50m isobath. This flow is partly captured in the northern part of the radar field (Fig. 2a). During November, the wind changes to northeasterly with larger speed. The surface currents demonstrate a typical winter pattern with a gulf-scale cyclonic circulation, an intense coastal current in the western part (Fig. 2 d,e), and several energetic jets with mean velocity of 0.3 m s -1 . The lack of river discharge in the model leads to decrease surface current velocity (by roughly 0.1 m s -1 ) in the vicinity of the Red river mouth in August (Fig. 2b,c) and large reduction of the southward coastal jet in November (Fig. 2f). the open-sea that occupies the southern area in August (Fig. 3b) to abundant thin filaments and small vortices across the gulf of the size from few to tens of km and with magnitude as large as (1 1.5)f (Fig. 3a,c). When the Rossby number Ro = ζ/𝑓 ~ O( 1), a departure from the geostrophic equilibrium in the flow is expected and smaller-scale processes can come into play, such as frontogenesis, mixed-layer instability [START_REF] Zhong | Submesoscale impacts on horizontal and vertical transport in the Gulf of Mexico[END_REF]. Consequently, these sub-mesoscale structures of the velocity fields can directly impact the dispersion process by increasing the rate of dispersion by an order of magnitude compared to the geostrophic turbulence [START_REF] Berta | Submesoscale Kinematic Properties in Summer and Winter Surface Flows in the Northern Gulf of Mexico[END_REF]. A significant difference in velocity gradients between the full model and the model without river discharge is the lack of a number of structures with high vorticity which are observed better in the vicinity of the Red river mouth in August and in November (Fig. 3b,d). Due to the coarse resolution of the HFR measurements, the vorticity field is much smoother, its magnitude is drastically reduced, and small-scale coherent structures are almost vanishing (results not shown).

The statistical characterization of the velocity field variability can be obtained from the surface turbulent kinetic energy spectrum, E(k) (Fig. 4). The spectra indicate a wide range of motions: from geostrophic and wind driven circulation at large scale to sub-mesoscale turbulent motions. The energy distribution appears different between two periods and higher for the GoT-Ref model configuration than for the GoT-noriv case.

In both models, the kinetic energy is concentrated at large scales (roughly at the length scale about 60 km) and decays more rapidly (steeper spectral slope ∼ 𝑘 -3 ) at scales below 20 km. In the range roughly above 20 km, despite some noise, the energy distribution follows 𝑘 -/3 law in both seasons. While the difference in energy distribution between two seasons is rather small, a significant difference is observed for model simulations with and without rivers forcing. In the latter case, the spectral slope from GoT-noriv model in summer appears steeper in a wide range of scales, 8 70 km, indicating a nonlocal dispersion regime from a Lagrangian perspective (Fig. 4a) while in November, above the range 10 30 km, the dispersion regime seems foster the local dispersion (Fig. 4b). tracking period, the surface currents are driven by the southwest monsoon while during the second period, the winter cyclonic circulation is formed, and the coastal current is reinforced as described in the previous section.

For particles released in August, trajectories are globally controlled by the cyclonic eddy located at the latitude 19 N which trapped the particles during the first 10 da s. The effect of particle trapping is clearl seen in trajectories drawn in the model and radar velocity fields (Fig. 5a,b). Ten days after the release when the eddy is weakened, particles separate in different directions in the model simulations. One group of them remains within the coastal zone close to XUAN radar site while another group is transported by the coastal buoyant jet southward, down to the gulf entrance (Fig. 5a). The direction of particle advection in the radar-derived velocity field is also southward for days 10 to 20, then it changes to eastward and northward (Fig. 5b). Agreement between the two transport patterns is found only for the first 15 to 20 days of tracking.

In November, a much better agreement in advection pathways in both HFR-derived and model velocity fields is found. In both cases, the trajectories of particles are related to the cyclonic circulation responsible for their transport southward. At the early stages, most of the particles are advected southward by the coastal jet. Later in time (roughly 10 days after the release), the particle drift direction changed to offshore, then reversed to northward after 15 days before continuing to the south after day 20. Regarding the number of particles remaining within the radar coverage area, 30% of particles left the area covered by HF radar measurements while this number is 10% for particles advected by the model velocity fields. The final location of the particles is shown in Fig. 5 in red color. At the end of the 30 days period, the particles advected in the radar velocity field show a more uniform distribution, while in the model simulation, the particles are separated into a number of different groups (Fig. 5). It is possible that larger dispersion in the model simulation is caused by irregularities in velocity fields due to sub-mesoscale motions which are probably resolved in the highresolution model.

To quantify the dispersion regime, we used relative dispersion and the FSLE. In the subsequent analysis of the turbulent dispersion, we used the estimation of the Rossby radius Rd performed by Nguyen-Duy et al. (2021).

In summer (April to September), large river discharge and thermal exchange with the atmosphere lead to the increase of stratification and Rd ranges from 10 to 15 km. During winter (Oct to Mar), strong wind and waves increase mixing in the water column. As a consequence, Rd has values of about 4 to 8 km and is even smaller in the shallow water area (see Fig. 4 in Nguyen-Duy et al., 2021). Therefore, in this study, the value of 8 km in November and 15 km in August are adopted for Rd.

Despite the differences in advection pathways, the relative dispersion regimes appear rather similar (Fig. 6).

For the cluster released in August, for both radar and model fields, relative dispersions indicate an exponential regime (exp (0.8t)) occurring during the first 10 days. It should be noted that in the case of the radar field, the pair separation grows slowly during the first 3 days, which can be explained by the coarser resolution of the HFR measurements (~6 km) that is larger than the initial separation (~2 km). During the next few days, the dispersion rapidly increases until reaching the distance equal to the grid resolution. The similarity in the exponential growth shown in Fig. 6a for both the HFR and the model indicates that relative dispersion during this initial period is dominated by processes with length-scale larger than the pair separation distance. On day 15, both curves flatten, suggesting that the particle pairs possibly reached eddies with typical length scale of about (20-30) km and larger (at later times). The evolution of relative dispersion follows a power-law behavior not far from Richardson regime for the model ( t 3 ) and ballistic regime for the radar ( t 2 ). The latter result is affected by larger statistical errors (Fig. 6a).

In November, the relative dispersion shows larger differences (Fig. 6b). Although the particles were released at the coastal current boundary which can rapidly separate the particles pairs, only the model indicates a high rate of exponential growth (exp(t)). This exponential phase is short, and last approximately 4 days. Compared to the model, the exponential growth in the HFR measurement is lower. After 10 to 15 days, both curves display a separation growth that is intermediate between a ballistic 𝑡 2 and Richardson 𝑡 3 regime, although it is more difficult to characterize the growth law for the HFR-derived data. Another approach used for identification of dispersion regimes involves the FSLE analysis (Fig. 6c,d). The FSLE uses the pair separation distance as an independent variable, thus, it is less affected by the superposition of contribution from pairs evolving under different dispersion regimes at the same time. In agreement with relative dispersion, the FSLE curves demonstrate a non-local dispersion regime for the scale range below 10 km, where the FSLE is constant, 𝜆 𝛿 ≃ 𝜆 , with 𝜆 the maximum Lagrangian Lyapunov exponent (Fig. 6c,d). In this range, the associated e-folding time 𝑇 1/ 2𝜆 is close to the values found in Fig. 6a, b, i.e., from 1 to 1.25 days in the model and from 1.25 1.6 days in the HFR velocity fields. The increase in exponential growth rate in correspondence with increased grid resolution (model vs HFR data) is consistent with previous works [START_REF] Poje | Resolution dependent relative dispersion statistics in a hierarchy of ocean models[END_REF] as pair dispersion at small separation is affected by small-scale motions resolved by the model.

At larger scales, somewhat common behavior can be recognized. In August, a transition of the dispersion regime is observed at the scale of about 30 km. Above this scale, the FSLE indicates a regime of dispersion close to Richardson one (λ(δ) ~ δ -2/3 ) for both velocity fields. A similar change of dispersion behavior is observed in November, but now at a scale of roughly 10 km, beyond which the FSLE is somehow steeper than Richardson regime and approaches the scaling 𝜆 𝛿 ∼ 𝛿 -1 typical of a ballistic regime.

Transport and dispersion of particles in the GoT from model simulations

Beside the coarse resolution of the HFR-derived velocity field, another weakness of the HFR measurements is the limited coverage which prevents the assessment of dispersion properties in an area in which the dynamics are strongly dominated by large-scale circulation features (gulf-scale eddy, coastal jet, etc.). Thus, to extend the study area to the whole region, we resort to numerical simulations. In particular, we are interested in assessing the effect of the river discharge which can rapidly generate sub-mesoscale structures and largely affect the dispersion in the small-scale range. A series of Lagrangian tracking experiments were performed during the same two periods (August and November) with particles seeded uniformly in the whole domain of model simulation. The average density of particle distribution at the release is 1.25 particles/km 2 . Two model configurations were used to assess the Lagrangian transport: the model with all forcings including the river discharge (GoT-Ref) and without river discharge (GoT-noriv). Fig. 7 shows the spatial distribution of particles 30 days after the release on August 3 (a, c) and November 10 (c, d) and their concentration per km 2 . In all cases, the spatial distribution of particles is found to be controlled by the flow field forced by the monsoon system. The particle distribution reveals a large difference between seasons. In August, the synthetic drifters are scattered throughout the gulf with a high density of particles observed in the middle and the northern gulf, while in November, most of the particles are located in the western gulf (Fig. 7b,d). At the end of the 30-day period, roughly 75% and 90% of particles released in August and November, respectively, remain inside the gulf. This suggests that the hydrodynamics of the gulf is quite local and the exchange with the open sea is limited. The number of particles leaving the domain is likely impacted by the open-boundary conditions. This effect will be evaluated in further studies.

Regarding the impact of the river discharge, the model simulation revealed that the distribution in August is quite similar in the southern gulf, whose dynamics are dominated by the inflow through the southern boundary (Fig. 7a,c). A much larger difference can be seen in the northern gulf. In the GoT-Ref case, the particles are organized into long filaments and wrapped around different eddies. Two regions of high particle density of particles are found there. The shape of the particle spatial distribution mimics the eddy-like structures with sizes that can reach 100 km. In the GoT-noriv case, the particle distribution is characterized by more disorder and homogenous scattering across the northern gulf. In general, the particle distribution in August reveals that mesoscale motions are largely enhanced by the river discharge which increase the number of well-organized structures characterized by higher concentrations of particles.

In November, both models show a large number of elongated structures in particle distribution field. These structures are associated with the unstable manifolds (attracting the materials drifting with currents) and thus structuring the transport of particles. They also indicate a cyclonic pattern in the GoT winter circulation (Fig. 7b,d). Moreover, the particle distribution in both model simulations shows higher concentrations along the middle Vietnamese coast, between 18° 20° N. In the GoT-Ref case, the lines of high-density particle distribution extend toward the south, demonstrating a tendency of particle alignment along the outer edge of the coastal current that carries them southward. In addition to the rapid southward transport, the fluctuation of the jet can facilitate the formation of eddies at the jet boundary, as described via idealized model [START_REF] Poje | Resolution dependent relative dispersion statistics in a hierarchy of ocean models[END_REF] which contribute to the enhancement of particle spreading.

In the GoT-noriv case, the majority of particles are concentrated along a line located roughly 50 km offshore in the northern part and in the southwestern part of the gulf, with large density (~5 particles per km 2 ) observed in the vicinity of the shoreline between latitudes 18°N and 19°N (Fig. 7d). 𝛿 2 ∼ 𝑡 3 at scales larger than 40 km (Fig. 8a).

In November, exponential growth in the GoT-Ref case with a similar rate occurs for a shorter timescale of 3 days (Fig. 8b) then quickly changes to a regime closer to the ballistic 𝛿 2 ∼ 𝑡 2 than to Richardson 𝛿 2 ∼ 𝑡 3 , hinting to particles separating due to the coastal buoyant jet. The ballistic regime holds until the end of the tracking period, when the largest separation of roughly 30 km is reached. The estimated slope for this period is 2.2. In the GoT-noriv case, exponential growth is observed is observed in a short early period, but with a smaller rate (~exp(0.6t)). The general features of the evolution of dispersion are not very different from those of the previous cases. The ballistic regime occurs after 5 days of drift. A possible explanation for the change to the ballistic regime in both configurations seems to be related to the particle entrainment inside the coastal current where the velocity shear causes a rapid separation in the cross-jet direction. In the whole, the rate of dispersion in the GoT-Ref case appears larger than that in the GoT-noriv simulation. The difference is caused by the presence of fine motions generated by the river plume advection where particles are stirred more effectively in the area that is not controlled by the mesoscale motions. The exponential growth of separation distance below Rd gives an e-folding time of 1.2 1.3 days for the GoT-Ref model and 1.7 2.2 days for the GoT-noriv configuration (Fig. 10), in fair agreement with the corresponding estimate obtained from relative dispersion (Fig. 8).

At larger scales (above Rd), a scaling slightly less steep than 𝛿 -2/3 and close to 𝛿 -1/2 is observed in August, likely due to relatively small-scale coherent structures that steepen the kinetic energy spectrum to 𝑘 -2 (see, e.g., [START_REF] Capet | Mesoscale to submesoscale transition in the California current system. Part I: flow structure, eddy flux, and observational tests[END_REF]. Note that in the GoT-noriv case the corresponding curve is less steep. Therefore, these results point to the enhancement of the particle spreading at scales comparable or smaller than 𝑅 , in the simulation including the river discharge. The largest difference between the two configurations is observed at small scales. It appears that including the river contributions enhances stirring and shortens the dispersion time by the factor of two in this range of scales (Fig. 10a).

The FSLE curves in November show rather different behavior. Above 20 km, the FSLE in both GoT-Ref and

GoT-noriv cases indicates a scale dependent dispersion rate in a broad range of scales. In GoT-Ref case, the FSLE scales as -1 (ballistic regime) until 70 km, which reveals a strong shear dispersion. This can be explained by the tendency of particles to be concentrated in the western gulf where the coastal current contributes to the rapid southward transport. At scales above 70 km, for both the full model and the GoT-noriv one, the FSLE is perhaps closer to -2/3 . In the larger scale ranges, the dispersion rate in November is on average slightly lower than in August (Fig. 10b). ). The LCSs shaped by these lines reveal a frontal zone of 25 km wide along the Vietnamese shore in both periods. Long and highly tangled ridge lines evidence the variability of the flow field there (Fig. 11a). Another region showing a large number of FSLE ridge lines is located further (~100 km) offshore. Most of the ridge lines in this region seem to be aligned along 50 m isobath (Fig. 11b).

In contrast to the full model, the LCSs derived from the GoT-noriv model simulation clearly lack many smallscale features in the coastal region. In addition, they indicate much weaker transport barriers (convergence zones) in the offshore region. FSLE ridge lines appear sparse and characterized by lower values compared to the full model simulation at the same period (Fig. 11c,d). Carrasco et al., 2012). We assessed the spatial variability of dispersion in the GoT by time averaging thirty FSLE maps for two periods. Several regions of high density of ridge lines can be seen in both seasons with the values of λ ranging from (0.2 to 0.6 days -1 ). Here we focus on four geographical regions: the northern gulf (A), the southern gulf (B), nearshore (C) and middle gulf region (D). It can be seen that regions of intense dispersion (λ = 0.6 days -1 ) are found in the northern gulf (location A) and along the shore (location B) characterized by also by the largest number of filaments and eddies (Fig. 12).

In August, thin ridge lines also indicate an intense dispersion in the north and middle regions (location A and C). Another interesting feature is a quasi-permanent ridge line with λ = 0.6 days -1 (location B) matching the location of a loop in the flow incoming from the southern boundary (Fig. 2b).

In November, the FSLE maps reveal two regions of intense dispersion located along the middle of the Vietnam coast (region C) and further offshore (region D) with λ ranging from 0.5 to 0.6 days -1 . The ridge lines in these two regions indicate the effect of the coastal jet and the branch of the current, originating from Hainan strait (Fig. 2e), on horizontal dispersion. Moreover, spreading of the ridge lines in location D (Fig. 12b) reveals large spatial variability of dispersion rate within an offshore band of 25-50 km large. LCSs can persist for up to several days thus limiting the exchange between the coastal and offshore regions located west and east of the ridge lines concentrated in sector D (Fig. 12b).

The intensity of stirring, which can be identified by taking the spatial average of the FSLEs over the whole area (d Ovidio et al., 2004), is appro imatel 0.4 da s -1 in the GoT-ref model. In the absence of the river discharge, the intensity of stirring decreases by roughly two times on average, down to 0.2 days -1 . These average values are reasonably in agreement with the FSLE plateau values λL at small scales from Fig. 10.

Discussion

In this study, we focus on the analysis of turbulent dispersion regimes and on the identification of LCSs in the GoT from realistic simulations using an approach based on particle pair dispersion. The analysis of relative dispersion is of great importance for the assessment of turbulent stirring in the open and coastal ocean at different spatial and temporal scales. Thus, the main objective of this study was to quantify the turbulent dispersion in the GoT, which is characterized by complex and multiscale dynamics, modulated by the monsoon wind and large buoyant input from a number of rivers (the largest being from the Red River). To the best of our knowledge, this problem has not been addressed yet in this region. Recent advances in the realism of the coastal modeling and observations using coastal radars allow the investigation of the turbulent dispersion using synthetic drifters advected in the radar and model velocity fields. Analysis of turbulent dispersion is performed for two periods, August and November 2015, representing typical summer and winter weather and river discharge conditions. Compared to other velocity data such as geostrophic currents from AVISO products, HFR measurements provide a near-real time, higher space-time resolution (6 km resolution, hourly data compared to roughly 27 km spatial resolution, daily data from AVISO). These data are more suitable for the quantification of the dispersion processes in the GoT. Our result indicates that relative dispersion statistics derived from the HFR tracking of particles are comparable to those obtained from the model. However, the spatial resolution of the HFR-derived velocity fields is still low compared to that of the model. As the short time and small-scale dispersion are sensitive to the details of the flow field, the coarse resolution of the HFR data makes the exploration of small-scale features difficult. In the small-scale ranges, the dispersion rate estimated from the radar measurements is significantly lower than the one from the model simulations in November (e-folding time of approximately 2 and 1 days, respectively).

Analysis of relative dispersion, relative diffusivity and FSLE from the numerical model indicated an exponential regime at sub-mesoscale which is consistent with estimations of these quantities for other ocean basins [START_REF] Beron-Vera | Statistics of simulated and observed pair separations in the Gulf of Mexico[END_REF][START_REF] Berti | Lagrangian Drifter Dispersion in the Southwestern Atlantic Ocean[END_REF][START_REF] Essink | Can we detect submesoscale motions in drifter pair dispersion[END_REF]van Sebille et al., 2015). At this scale, the e-folding time estimated using the relative dispersion and the FSLE is similar for both seasons but appears strongly dependent on the river runoff (e-folding time of roughly 1.25 days for the full model and 2.5 days in the GoT-noriv model simulation). At larger scales, the dispersion exhibits a scale-dependent behavior indicative of either Richardson local dispersion or ballistic separation, controlled by strong mean currents.

Such a transition between regimes occurs at scales close to the Rossby radius of deformation. In general, the pair dispersion properties in the GoT, obtained from the analysis of three Lagrangian diagnostics (relative dispersion, relative diffusivity and FSLE), are consistent with the results of previous studies, based on the analysis of the global drifter dataset [START_REF] Corrado | General characteristics of relative dispersion in the ocean[END_REF].

It is argued that a super diffusive regime can occur in the presence of strong shear and is characterized by the power law behavior ~𝑡𝛾 𝛾 1 , as reported in some other studies performed in different regions: in the Adriatic Sea, is found 1.8 (Ha a et al., 2008), in the Gulf of California, 1.9 (Zavala Sans n, 2015). In our case, is found larger, 2.2. However, rather clear hints to Richardson regime 𝛾 3 are also found over some scale subranges for both model configurations, with and without river discharge.

In regions characterized by strong sub-mesoscale dynamics, different limiting FSLE values were found for scales δ < 10 km: λL, ~ 1 1.1 days -1 in the Gulf of Mexico, [START_REF] Choi | Submesoscale dynamics in the Northern Gulf of Mexico. Part III: Lagrangian implications[END_REF], or λL ~ 0.6 days -1 in the Liguro basin [START_REF] Schroeder | Relative dispersion in the Liguro-Provençal basin: From sub-mesoscale to mesoscale[END_REF]. In the Gulf of La Spezia, the value of λL can readily reach 2.5 days -1 [START_REF] Schroeder | Targeted Lagrangian sampling of submesoscale dispersion at a coastal frontal zone[END_REF]. In our case, the value of λL for the same scales (δ < 10 km) is clearly weaker, ranging from 0.3 to 0.4 days -1 . This indicates that the contribution of sub-mesoscale motions in turbulent dispersion in the Gulf of Tonkin is weak. This seems reasonable because the GoT is a quite shallow basin with limited potential energy which can be stored in the mixed layer where small eddies can be hashed by the wind. Our results are in agreement with those previous studies which reported and demonstrated that dispersion at the sub-mesoscale can be induced non-locally by mesoscale motions (e.g., Haza et al., 2008[START_REF] Poje | Resolution dependent relative dispersion statistics in a hierarchy of ocean models[END_REF]. However, it should be noted that the model resolution in the GoT (1 2 km in the offshore area to 4.5 km at the open ocean boundary) probably prohibits the exploration of the local stirring generated by small-scale turbulent motions.

The impact of the river discharge on the surface dynamics of the GoT is explored using the model simulation with and without rivers runoff. During both periods considered, the river input in the model enhances the dispersion in a wide range of scales, from sub-mesoscale up to 50 km. The river discharge causes a significant change of e-folding time, which decreases from roughly 2.5 days to 1.25 days. The river discharge (when introduced in the model) creates buoyancy gradients and enhances the sub-mesoscale dynamics which in turn modifies the Lagrangian statistics. However, the difference between two numerical experiments, with and without river discharge, was not found to significantly alter the dispersion scaling behaviors. A possible reason could be that we examined statistics averaged over all particle pairs in the whole domain, while the river discharge only enhances the dispersion close to the western coast and locally affects the particles released there.

Lastly, we quantified the horizontal stirring caused by small-scale motions using the spatial distribution of FSLEs in different seasons. The LCSs, identified from backward FSLEs using the model velocity field, were analyzed to assess the horizontal stirring by small-scale processes, and appear to be very sensitive to the river forcing. Moreover, the FSLE maps clearly show the regionality of small-scale motions and their impact on dispersion processes. In August, a large concentration of FSLE ridge lines is found in the northern gulf and in the southern gulf where the transport and dispersion are strongly influenced by the incoming flow from the southern open boundary. During August, when the river discharge is at its peak, the river plume is driven northward and can detach from the coast due to the persistent southwest wind (Nguyen-Duy et al., 2021). The surface layer is then decoupled from the subsurface layer. This creates favorable conditions for the development of sub-mesoscale motions via roll-up instabilities of larger scale structures (Klein et al., 2011;[START_REF] Roullet | Properties of steady geostrophic turbulence with isopycnal outcropping[END_REF][START_REF] Berti | Lagrangian pair dispersion in upper-ocean turbulence in the presence of mixed-layer instabilities[END_REF]. Alternatively, sub-mesoscale motions could be also generated by a complex topography in the northern gulf (small islands, complex bathymetry, etc.). This can explain a high concentration of large FSLE values in the vicinity of the Red river in the GoT-Ref model (location A in Fig. 11a). In winter, the spatial distribution of FSLEs reveals the existence of another region of high concentration of FSLE ridge lines matching the outer limit of the coastal current which strongly affects the dispersion processes by organizing the particle distribution.

Summary

Relative dispersion in the coastal ocean environment, characterized by a wide range of scales of flow variability is still largely unknown. At scales less than a few km, the dispersion processes have been experimentally investigated using surface drifters in a number of limited size coastal regions. The experimental assessment of turbulent dispersion at larger scales is difficult to realize due to the cost of such experiments. In the present study, we adopted an alternative approach and investigated the turbulent dispersion and transport pathways in the GoT by using a large number of virtual drifters advected in the velocity fields derived from a highresolution numerical model and HFR measurements.

For both types of fields, the results revealed a high degree of consistency regarding the dispersion regime at scales larger than 20 km. At smaller scales, the dispersion rate in the radar velocity field appeared two times lower which is explained by the comparatively low spatial resolution of radar measurements.

The Rossby radius of deformation, Rd, is a key lengthscale that separates the regimes of turbulent dispersion.

It was estimated from numerical simulations in the GoT performed by Nguyen-Duy et al. (2021). Our results revealed a non-local dispersion regime at scales below Rd and local dispersion at larger scales, which is in agreement with the results of earlier investigations performed different world oceans (e.g., Haza et al., 2008;[START_REF] Corrado | General characteristics of relative dispersion in the ocean[END_REF].

The seasonal variability in transport patterns and dispersion regimes was further characterized. The time evolution of dispersion showed an exponential growth during first 5 8 days of drift for both seasons, followed by a power law regime during the next 6 to 20 days. The ballistic regime was found to dominate the dispersion at mesoscales in winter while Richardson super diffusion regime was diagnosed in summer. This change in regime is mainly related to a combined effect of river discharge and summer monsoon wind which enlarges the river plume area and enhances stirring in the coastal region extending up to 100 km offshore. In November, the ballistic regime is probably induced by strong northeastern monsoon wind providing the coastal current strengthening and local shear dispersion intensification. The non-local dispersion regime (in sub-mesoscale range) was not clearly identified, suggesting that the mesoscale circulation in the GoT largely governs the dispersion at small scales.

Four regions with high level of dispersion were identified in the GoT: a nearshore region under the fresh-water influence, an offshore region distant by ~100 km from the coast, a central part of the northern gulf, and a region close to the southern gulf entrance. A seasonality in dispersion strength, caused by the variability of wind forcing and freshwater input by the Red river, was detected.

Summarizing, the Lagrangian diagnostics appeared very useful for characterization of turbulent dispersion regimes and identification of regions of high stirring. This is of primary importance for applications involving multi-scale interactions in the coastal ocean, the transport of particulate biological or anthropogenic material.

Figure 1 :

 1 Figure 1: Location of the Gulf of Tonkin. Radar sites providing the data used for the present analyses are indicated by triangles. Black dots indicate the HF radar coverage. The dot spacing is 12 km. Color shading represents a snapshot of the Chl-a concentrations (in log10 scale) obtained from Aqua/MODIS on September 24, 2015. Extension of high Chl-a concentration patterns from the nearshore blooming region seaward is visible through the elongated filaments and small eddies.

  Eq. (3), with 𝛿 ≃ 1.86 km and 𝛿 𝑓 13 km. Such values of 𝛿 and 𝛿 𝑓 allow to adequately visualize the relevant coherent structures and transport patterns associated with small-scale processes.

Figure 2 :

 2 Figure 2: 30-day mean surface currents in the GoT from HFR (a, d) and numerical model in two configurations: with river discharge (b, e) and without river discharge (c, f) in August and November 2015.

Figure 3 :

 3 Figure 3: Snapshots of the relative vorticity normalized by Coriolis frequency in August and November 2015 from model simulation with (a, c) and without river discharge (b, d). Gray line shows 50 m isobath.

Figure 4 :

 4 Figure 4: Power spectra of surface horizontal turbulent kinetic energy in August (a) and November (b) for two model configurations: with and without river input.

Figure 5 :

 5 Figure 5: Trajectories of the subset of particles released close to the XUAN radar site (black squares) during both periods of study from the model (left panel) and from the HFR data (right panel) for 30 days since release. The drift time is color-coded. The particle final positions correspond to the red points.

Figure 6 :

 6 Figure 6: Scaling properties from relative dispersion (a, b) and FSLE (c, d) for particles advected by model and HFR-derived velocity field in August and November. The horizontal (a, b) and vertical (c, d) lines correspond to separation distances equal to the Rossby radius 𝑅 .

Figure 7 :

 7 Figure 7: Particle distribution in August and November 30 days after the release from the full model (GoT-Ref) (a, b) and without river discharge (GoT-noriv) (c, d) (unit: number of particles per km 2 )

Fig. 9

 9 Fig.9shows relative diffusivity κ as a function of the separation distance δ in August and November respectively. The variability of this quantity is more affected by noise as κ is the time-derivative of relative dispersion. Thus, following[START_REF] Koszalka | Relative dispersion in the Nordic Seas[END_REF], we binned κ(δ) in such a way that a given value of δ is a multiple of the previous one, as one would do for the computation of the FSLE. Only the values with more than 50 data points for each bin are considered in the analysis. Despite some wiggles, the scaling properties of κ are consistent with the relative dispersion ones. During August and November, κ scales as δ 2 in the range below Rd. This indicates the dominance of non-local dispersion at these scales. The transition from the nonlocal to the local regime (κ ~ δ 4/3 ) is observed for scales close to Rd (and in fact slightly smaller). The estimated exponent ranges from 1.2 to 1.5 (Fig.9), and it is smaller in November. The comparison of diffusivity demonstrates that the freshwater runoff increases the diffusivity in the whole domain.

Figure 8 :

 8 Figure 8: Relative dispersion as a function of time from the release in August (a) and November (b) in linearlog scales. In this case, particles were released uniformly in the GoT domain. The horizontal lines correspond to separation distances equal to the Rossby radius Rd.

Figure 9 :

 9 Figure 9: Relative diffusivity as a function of the separation distance for the two tracking periods in August (a) and November (b) and for the two model configurations (with and without river discharge). In this case, particles were released uniformly in the GoT domain. The vertical lines correspond to separation distances equal to the Rossby radius 𝑅 .

  Figure 10: FSLE λ(δ) versus the separation distance in August (a) and November 2015 (b), for particles released uniformly in the GoT domain. The vertical lines correspond to separation distances equal to the Rossby radius 𝑅 .

  Figure 11: FSLE distribution for two particular dates on September 02 and December 10 using two model simulations. Gray line shows 50 m isobath.

Figure 12 :

 12 Figure 12: 30-day mean FSLEs for the period of August and November from the two model simulations: GoT-Ref (a, b) andGoT-noriv (c, d)

Acknowledgments

This stud was supported b the Universit of Littoral C te d Opale and the R gion Hauts-de-France PhD fellowship. Chl-a data were provided by ACRI, GlobColour products (https://hermes.acri.fr/). We would like to thank Tran Trung Kien for help with the data. A ho con ib ion M.C.T, A.S, S.B contributed to the conception and scientific design of the study. M.C.T performed the data analysis and wrote the original version of the manuscript. A.S, S.B and N.A deeply involved in the analysis,

Funding

This study was funded by Université du Littoral-C te d Opale for PhD scholarship fellowship.

manuscript organization, and co-writing. N.A and T.ND developed and provided the numerical model results.

All authors actively contributed, comments and revised the manuscript.

Conflict of interest

The authors declare that they have no conflict of interest.


Vinicius Beltram Tergolina

Enrico Calzavarini

Gilmar Mompean

Stefano Berti

Surface light modulation by sea ice and phytoplankton survival in a convective flow model

Plankton dynamics depend in a complex manner on a variety of physical phenomena, according to both experimental and numerical data. In particular, experimental field studies have highlighted the relation between phytoplankton survival and turbulent upwelling and downwelling from thermal convection. Recent numerical works have also shown the importance of accounting for advective transport by persistent structures in simulation models. In nutrient-rich polar marine environments phytoplankton blooms are critically limited by light availability under ice-covered waters. Such heterogeneity of the light intensity distribution, in association with a large-scale coherent fluid flow, can give rise to nontrivial growth dynamics. In this work we extend a previous advection-reaction-diffusion model of phytoplankton lightlimited vertical dynamics in the presence of convective transport. Specifically, we consider horizontally heterogeneous light conditions through the use of two regions with different production regimes, modelling the absence (presence) of light under (in between) obstacles. Such a model is intended as an idealized representation of nonuniformly ice-covered polar waters. By means of numerical simulations, we find that the main role of advective transport is to hinder phytoplankton growth, but also that such effect depends on the positions of the obstacles with respect to the upwelling and downwelling flow regions. Furthermore, we show that the sinking speed due to the density difference between phytoplankton organisms and water, while small, plays an important role, which depends on how it adds to the flow. These results indicate that advective transport can have a crucial impact on the survival conditions of sinking phytoplankton species in polar environments.

Introduction

In phytoplankton ecology, biological and fluid dynamics are closely linked, ruling how planktonic populations interact with their environment and other organisms. Interdisciplinary work aimed at studying the interplay between plankton biology and the concurrent physical factors has allowed to considerably refine the understanding of the subject [START_REF] Mann | Dynamics of marine ecosystems: biological-physical interactions in the oceans[END_REF][START_REF] Kiørboe | A mechanistic approach to plankton ecology[END_REF][START_REF] Prairie | [END_REF]. Notably, both modeling studies (see, e.g., [4][5][6]) and direct measurements [7] demonstrated the importance of vertical advective transport, underlining the relation between turbulent convective fluid motions and phytoplankton growth dynamics.

Polar marine ecosystems are the seat of important plankton blooms (see [8] for a review). Understanding the dynamics of the latter is key to correctly estimate primary productivity [7,[9][10][11], particularly under a changing climate scenario. In the Arctic ocean, several studies reported under-ice phytoplankton blooms that were initiated by the onset of ice melt [7,9,11,12]. Nevertheless, it is still not known whether under-ice blooms are controlled by the increase of light transmittance, leads (openings in the ice), upwelling from springtime convective mixing, biological factors, or other factors (most likely a combination of more than one). Interestingly, however, it was observed that waters beneath sea ice with leads had significantly lower phytoplankton biomass, despite high nutrient availability, than regions without (or with few) leads [7]. It was argued that this phenomenon is controlled by the onset of convection below leads, which enhances mixing to larger depths, and reduces the travel time along the vertical. Consistently with Sverdrup's arguments [13], then, the bloom should be suppressed due to the mixed-layer depth being larger than the critical depth (where the integrated production is balanced by the integrated loss). Note that convection is here to be intended both as a thermo-solutal buoyancy-induced phenomenon, due to brine rejection from the ice and surface heat losses, and as wind-induced mixing in leads [7]. It is worth remarking, too, that phytoplankton blooms beneath loosely consolidated ice will likely become more common in the future Arctic ocean, as lead formation becomes more frequent due to increasingly thinner and dynamic ice coverage [14][15][16]. This could significantly alter primary production, and have important consequences on local marine food webs.

In this work we take an approach based on idealized models, to explore the effect on phytoplankton dynamics of the horizontal modulation of light intensity at the surface of a vertical fluid layer. For this purpose, we extend a previous two-dimensional (2D) advection-reaction-diffusion model that we used to investigate the role of large and small scale fluid motions on the same biological system, but with uniform surface light intensity [6]. The role of vertical turbulent diffusivity on light-limited phytoplankton dynamics has been extensively studied in the past [17][18][19][20] in onedimensional (1D) systems. Our model can be seen as a 2D generalization of the one introduced in [19]. Moreover, we also consider the presence of a coherent large-scale flow (as in [6]) and of horizontally varying incident light intensity at the surface. The flow field will be given in terms of a streamfunction corresponding to a kinematic cellular flow. While this is clearly a strong simplification with respect to a truly dynamic velocity field, it allows us to use realistic biological parameters and domain size. The cellular flow structure is intended as a simplified description of buoyancy and wind-driven flows [START_REF] Thorpe | The turbulent ocean[END_REF], such as convective motions and Langmuir circulations [START_REF] Craik | [END_REF][23][24], often encountered in the upper ocean. Beyond the environmental implications discussed above, from a more fundamental point of view, based on the overwhelming importance of the large-scale flow and of its downwelling region that we found in a previous work [6], in this study we wanted to see how phyotplankton dynamics change when reproduction is forbidden over either the upwelling or the downwelling flow region. As also in this case we expect a weak effect of small-scale flow features, we only consider a steady large-scale (cellular) flow. We then place a light-blocking obstacle, either over the downwelling region or over the upwelling one.

Through this approach we aim at assessing the importance of the different transport mechanisms on phytoplankton survival, in the presence of such vertically and horizontally heterogeneous growth dynamics. We also examine different values of the sinking speed, due to phytoplankton being typically heavier than water. This corresponds to considering different algal species, and appears interesting on the basis of the experimental evidences about the varied species composition of phytoplankton in under-ice blooms [9].

This article is organized as follows. The model is introduced in Sec. 2. In Sec. 3 we report the results of our numerical study, and we discuss the impact of the different transport terms, by separately considering them in a systematic way. The main outcome and the conclusions are presented in Sec. 4.

Model

Our model is intended as an idealized representation of phytoplankton dynamics in waters that are partially covered by ice. Specifically, we consider light-limited growth in a 2D vertical fluid layer, in the presence of obstacles at the surface mimicking ice blocks that do not allow the transmission of the incoming radiation to the water underneath. A schematic is provided in Fig. 1. While our description is an over-simplification of a realistic under-ice ecosystem, the present setup can provide indications useful to develop more complex biophysical models, and is in principle quite easily reproducible in laboratory experiments.

To describe the dynamics of the system is then necessary to specify the evolution equation of the population density field θ(x, z, t) (number of individuals per unit volume) in both the ice-free regions and in the ice-covered ones. In the regions receiving light (where phytoplankton can reproduce), we have:

∂θ ∂t = [p(I) -l] θ -v • ∇θ + D∇ 2 θ, (1)
as in [6], where p(I) is the light-dependent production rate and l the loss rate. In the regions below obstacles, instead, due to the absence of light, p(I) = 0 and one has:

∂θ ∂t = -lθ -v • ∇θ + D∇ 2 θ. (2)
The fluid layer in which we examine these dynamics is intended to represent the mixed layer (the water layer of nearly constant turbulent intensity found in upper oceans and lakes [START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF]). Its horizontal and vertical sizes will be denoted L x and L z , respectively. The top and bottom boundaries are assumed to be rigid (as in [6,19]).

The production term in Eq. (1) accounts for both water background turbidity, via the coefficient κ bg , and population self-shading, with an attenuation factor κ. Its functional form is:

p(I) = p max I H + I , (3)
where p max is the maximum specific production rate, H is a half-saturation constant and the time-and depthdependent light-intensity is expressed as follows, according to Lambert-Beer's law: with I in the incident light (at the surface, where z = 0). Note that we will assume I in = 0 in correspondence with the surface obstacles. The biological parameter values used in our study, representative of realistic situations [START_REF] Huisman | [END_REF], are reported in Table 1. They correspond to average measured values, with growth parameters obtained from freshwater phytoplankton species, and κ bg from clear lakes and coastal areas [20]. The coefficient D in Eqs. ( 1) and ( 2) represents an effective diffusivity, due to both small-scale unresolved turbulent motions and possible non-directed swimming behavior. Advection is realized by a 2D incompressible flow u = (u x , u z ) and phytoplankton is assumed to sink with a speed v sink ẑ, where ẑ is the unitary vector pointing downward in the vertical direction; the total velocity appearing in the evolution equations is thus v = u+v sink ẑ. For the phytoplankton sinking speed we take the value in Table 1 as the reference one, and we will later vary it.

I(z, t) = I in e -z 0 κθ(s,t)ds-κ bg z , (4)
For the velocity field we choose a cellular flow, which was originally introduced as a kinematic model of thermal convection fluid motions [27]. Its expression is obtained as u = (-∂ z Ψ, ∂ x Ψ ) from the following streamfunction:

Ψ (x, z) = - U k sin(kx) sin(kz), ( 5
)
where U is the maximum velocity and k = π/L z a wavenumber. The flow pattern then corresponds to two counterrotating steady vortices of size L z , in the (x, z) plane. It is worth stressing that this flow comprises both a (central) downwelling region and two (side) upwelling ones, which typically play a relevant role on phytoplankton dynamics [5,6]. Note that eddying motions on scales smaller than L z , and unsteadiness, could be in principle taken into account by adding to the streamfunction of Eq. ( 5) analogous terms with larger wavenumbers, and an explicit time dependency (see, e.g., [6,28,29]). Our choice to neglect them, here, is motivated by a previous study without the surface light modulation [6], where we found that small unsteady eddies play a secondary role. Indeed, in that case, the fate (survival or extinction) of phytoplankton was primarily determined by the large-scale coherent component of the flow.

In the following we will consider two different geometrical configurations, corresponding to light-blocking obstacles, of fixed total size, positioned either above the downwelling flow region (LBD case), as in Fig. 1, or above the upwelling one (LBU case). The horizontal and vertical domain sizes will be respectively set to L x = 60 m and L z = 30 m.

Phytoplankton sinking velocity

Considering the chosen biological parameters (Table 1), a mixed layer depth of 30 m is smaller than the critical depth (about 50 m). Hence, in the absence of surface light heterogeneity and of a fluid flow, with the values of diffusivity explored, 1 cm 2 s -1 ≤ D ≤ 20 cm 2 s -1 , which appear quite realistic for the ocean [6,19,[START_REF] Huisman | [END_REF], a bloom occurs. The presence of obstacles considerably reduces production, due to the absence of light below them. However, we made sure that also in this case phytoplankton blooms when U = 0. This allowed us to investigate the impact of advection on the dynamics of the system, including a possible transition from population survival to extinction. Remark, too, that working with a relatively small system, but still apt to capture the salient dynamics, is computationally convenient. We numerically study the dynamics of Eqs. ( 1) and ( 2), which are simultaneously integrated in our rectangular domain (with L x = 2L z ) by means of a pseudo-Lagrangian algorithm [30,31] (see [6] for more details). We adopt periodic and no-flux boundary conditions along the horizontal (x) and the vertical (z), respectively (as in [6]). The initial condition is a low uniform population density [θ(t = 0) = 5.5•10 6 cells m -3 ]. To analyze the blooming conditions of the phytoplankton population we rely on the temporal behavior of the average biomass density,

θ (t) = 1 L x L z Lx 0 Lz 0 θ(x, z, t) dxdz, (6)
and the per-capita growth rate,

r p (t) = 1 θ ∂ θ ∂t . (7)
The per-capita growth rate is expected to attain a statistically constant value r p (after an initial transient), corresponding to exponential growth (r p > 0) or decay (r p < 0) in the early regime before the onset of nonlinear dynamical effects.

Results

We begin our numerical study by exploring the dynamics in two partially illuminated setups, both in the absence and in the presence of the fluid flow. As explained in Sec. 2, the blocking of light at the surface is realized by introducing obstacles of zero light transmittance, with different positions according to the cases (LBD, LBU) considered. Given the large number of control parameters (for both biology and transport) in the model, we decided to work with a single light-transmittance value associated with the obstacle. We then chose to use fully absorbing obstacles, in order to increase the contrast between the zones where light availability permits reproduction and those where this does not occur. Such a simplifying choice is also intended to ease the interpretation of the effects of the different transport and production regimes. Indeed, as we shall see, already in the present case the dynamics are quite complex. The light transmittance of ice is nevertheless variable and is a subject of debate (see [32] for recent measurements of ice thinning in the Arctic ocean). Some snapshots of the field θ at a given time are shown in Fig. 2, for both the LBD and LBU cases, and for increasing values of the flow intensity U . This figure is obtained with the reference sinking speed v sink = 0.04 m h -1 [6,19]. When advection is not present the two light-modulation cases give identical results, as expected. With nonzero advection, instead, important differences appear. While in the LBD case phytoplankton accumulates closer to the surface and presents small differences with respect to increasing flow intensity, in the LBU case the population spreads more in the vertical, reaching larger depths. By measuring the per-capita growth rate r p versus time for different values of U , we could assess that phytoplankton growth is hindered by advection, independently of the light-modulation considered. Indeed, the (constant) value of r p reached at large enough times, monotonically decreases with increasing U in both the LBD and LBU cases [see Fig. 3, where r p is normalized by the total (birth minus death) growth rate at the surface

r b = p max I in /(H + I in ) -l].
This effect parallels the one we already found in the absence of obstacles [6]. However, the decrease of r p in the two situations is quantitatively different, and indeed larger when light does not penetrate in the upwelling region, so that an LBD/LBU asymmetry emerges. In the remainder of this work, we focus on the understanding of this feature.

We can first observe that the orientation, with respect to the light gradient, of the fluid velocity alone (which is opposite in the LBD and LBU cases) cannot be responsible for this effect. To see this, one can rewrite Eq. (1) as follows:

∂θ ∂t = g(x)θ -u • ∇θ, (8)
where g(x) indicates the spatially dependent growth rate, and diffusion and sinking have been neglected. Note that in our system the growth rate depends on both the vertical and the horizontal coordinates, due to the external light blocking and to self-shading. The fluid velocity u(x) depends on space, too, but not on time. In a Lagrangian reference frame the above equation reads:

θ = g(x(t))θ, ẋ(t) = u(x(t)). (9)
Given an intial value for θ at time t = 0 and a final position x(t), the solution of Eq. ( 9) is:

θ(x(t), t) = θ(x(0), 0) e t 0 g(x(t ))dt (10) = θ(x(0), 0) e |x(t)| |x(0)| g(x ) |u(x )| dx . ( 11
)
It is therefore clear that the growth of θ along a given Lagrangian trajectory (coinciding with a streamline, due to flow stationarity) depends on the velocity at which the path is travelled but not on the orientation of the velocity itself. Moreover along a closed path, like a cellular-flow streamline, the variation of the reactive-scalar density will be independent of the circulation direction. A full loop starting with downward flow and followed by an upward flow will produce the same phytoplankton growth if it occurs in the reversed order. The presence of a sinking speed, which modifies the Lagrangian equation of motion to ẋ(t) = u(x(t)) + v sink ẑ, can change the scenario because a descending flow path will be covered faster than an ascending one, leading to a smaller growth for the first than for the latter. Note that the presence of sinking also has an effect on the shapes of planktonic trajectories, in particular not all the trajectories in the cellular flow will be closed [33]. It is therefore the (slow) sinking movement that is responsible for breaking the symmetry of the growth dynamics between the upwelling and downwelling flow regions. The presence of diffusion will also act in a similar manner by allowing parcels of population density to cross the flow streamlines.

Figure 4 reports the measure of the per-capita growth rate -normalized by the total (birth minus death) growth rate at the surface r b -versus time, for different values of the sinking speed and a fixed flow intensity. Panel (a) in the figure refers to a small value of the diffusion coefficient, while panel (b) is for a larger one. The role of diffusivity is quite subtle and we will come back to it later in this section. Independently of this, the results clearly show that in the absence of sinking (v sink = 0) the system is insensitive to the position of the obstacle. The value of r p is in fact the same, within numerical accuracy, in the LBD and LBU cases. However, it decreases in the presence of sinking when light is blocked over the upwelling region, as it may be expected. It is perhaps more surprising that it instead increases in the LBD case (which is more difficult to see for large D in Fig. 4b due to the scale). For small diffusivity, it can even turn a no-bloom situation (r p < 0) into a bloom one (r p > 0). The combined effect of fluid advection, sinking and diffusion can be appreciated through the following argument, which helps rationalizing the previous observations. We consider the transformation (see also [34]):

θ(x, t) = exp w • x 2D θ(x, t) (12)
where w is a, yet to be defined, constant (in both time and space) velocity. Plugging Eq. ( 12) into Eq. ( 1) we get:

∂ θ ∂t = p(I) -l - w • v 2D + w 2 4D θ -(v -w) • ∇ θ + D∇ 2 θ. ( 13
)
Taking w = v sink ẑ, and using the fact that v = u + v sink ẑ, we obtain

∂ θ ∂t = p(I) -l - v sink u z 2D - v 2 sink 4D θ -u • ∇ θ + D∇ 2 θ. ( 14
)
This is the evolution equation of a field advected by the fluid flow u, where sinking is not explicitly present, but with a modified reaction term that now depends on the sign of the vertical velocity component u z (and also on space).

It is thus clear that the local growth rate, i.e. g = p(I) -l -v sink u z /(2D) -v 2 sink /(4D), is always smaller in the downwelling regions (u z > 0) than in the upwelling ones (u z < 0). In short, g d < g u , where d, u label down-and upwelling regions, respectively. We also note that the quadratic term in v sink can be neglected nearly everywhere, because in practice v sink |u z | ∼ U (this may be not true close to the top and bottom boundaries, where u z also becomes very small). The above reasoning is independent of the presence of the obstacle blocking light at the surface. When the latter is present we can expect that the LBU situation is characterized by a growth rate similar to the one observed in downwelling regions, g LBU ∼ g d and, on the opposite, for the LBD case: g LBD ∼ g u . This implies that g LBD > g LBU , in agreement with the observed trends.

One may then be tempted to extrapolate this result to the per-capita growth rate, r p = ∂ t θ / θ , to obtain the following predictions for the two situations considered:

r p,LBD ∼ p(I) L -l + 2v sink U π 2 D - v 2 sink 4D , (15)
r p,LBU ∼ p(I) L -l - 2v sink U π 2 D - v 2 sink 4D , (16)
which, however, assume a weak spatial coupling between the growth rate g and θ, and where we used the average vertical velocity in the regions receiving light, u z L = ±(2/π) 2 U . More importantly, here we have neglected any other effects due to both advection and diffusion. In particular, this argument completely ignores the lateral transport, between the illuminated and shaded regions, due to both advection and diffusion. The flux associated with this process is u x θ -D∂ x θ, with u x = ±U cos(πz/L z ) at the boundaries between the two regions, and it should be included to improve the present prediction. Indeed, the expressions in Eqs. ( 15) and ( 16) fail to describe the behavior of the large-time value of r p as a function of v sink , unless some fitting coefficients are introduced to account for the missing physics. While for small diffusivity the fitting coefficients are comparable in the LBD and LBU cases, which may suggest that the missing information resides in the z-dependence of u z and in the appropriate choice of the region where to average, this is not the case for large D. In fact, in this case, a striking feature that is already apparent from Fig. 4b is the asymmetric response of the system with respect to the position of the obstacle. This fact cannot be explained by Eqs. ( 15) and ( 16), which are symmetric under the exchange of U into -U . Taking into account the role of lateral transport between shaded and illuminated regions analytically is not trivial. Remark, too, that in the large D limit analytical predictions are not available even in the 1D no-flow (fully illuminated) model [19], starting from which ours was constructed. In the following we then resort to an argument based on Lagrangian paths to provide a qualitative explanation of the asymmetric response of the system.

From a Lagrangian point of view, in the presence of diffusion, the complete equation of motion describing fluid transport is ẋ = u(x(t))+v sink ẑ+ √ 2Dη, where the white noise term √ 2Dη represents diffusion. Sample distributions, at an arbitrary time t > 0, of N = 2000 particles evolving according to this equation are shown in Fig. 5 for D = 1 cm 2 s -1 (red circular points) and D = 10 cm 2 s -1 (blue square points). The two panels correspond to the LBD (Fig. 5a) and the LBU (Fig. 5b) case. In each of them, the particles were initially located in a small area close to the top of the illuminated region (green points). As expected, for larger D particles spread more into the domain. It is here worth emphasizing the lateral spreading: particles originally in a well-lit region are much more likely to reach a dark region, where reproduction is not possible, if their diffusivity is larger. In the LBD case this implies that, when D is large enough, lateral diffusion can overcome the action of the sinking velocity, which helped maintaining particles in the upwelling reproductive region. In the LBU case, particles are also more easily brought to deadly regions by increased diffusivity, but this occurs while sinking contributes moving them downwards, where growth is also importantly reduced due to light absorption along the vertical. Therefore, the overall effect is weaker than in the previous case. This explains the asymmetry observed in the behavior of the per-capita growth rate as a function of the sinking speed, with respect to the position of the obstacle blocking light. Indeed, in the LBD case, the positive action of v sink for growth is only barely visible for D = 10 cm 2 s -1 in Fig. 4b, as it is masked by the impact of diffusivity. On the contrary, in the LBU case, such a large value of the diffusion coefficient does not change the general phenomenology (with respect to the D = 1 cm 2 s -1 case of Fig. 4a), as this was already unfavorable for growth, due to downward advection by both the sinking and vertical fluid velocities in the illuminated downwelling region. We conclude this discussion by noting that the same reasoning can be made also considering transport backwards in time, which allows interpreting the results in terms of the generalization of Eq. ( 10) accounting also for diffusive transport:

θ(x, t) = θ(x(0), 0) e t 0 g(x(t ))dt η . (17)
The above equation (see [30] for more details) states that the population density at time t > 0 and at position x is constructed by advancing in time the reactive dynamics, starting from the value of the field θ(x(0), 0) at the Lagrangian origin x(0) of a trajectory ending at x at time t, and averaging over noise (η) realizations. Clearly, when the diffusion coefficient increases many more Lagrangian paths reaching x at time t originate from regions where light cannot penetrate, bringing already dead material at the considered position at a later time. In the LBD case, the effect is considerable, as the spatially and temporally varying per-capita growth rate g(x(t)) in Eq. ( 17), at a location x in the illuminated (and upwelling) region, is now determined by the many trajectories coming from dark regions, which was not the case for smaller D. In the LBU case, instead, also with small values of D the contributions to g(x(t)), with x in the illuminated (and downwelling) region, are essentially associated with trajectories starting in regions where light is blocked and the population has already died. Therefore, increasing D does not alter the qualitative behavior of the system in this case.

Conclusions

We numerically investigated the dynamics of sinking phytoplankton in a stirred fluid layer in which light availability is spatially inhomogeneous, along both the vertical and the horizontal. Our system is assumed to be illuminated from above, as is typical in aquatic ecosystems receiving energy from the incoming solar radiation. Vertically, light intensity decreases due to absorption by water and phytoplankton self-shading. Horizontally, its penetration from the surface into the interior of the medium is critically affected by the presence of opaque obstacles at the top boundary. We focused on the impact of advection, and more generally of the different transport processes occurring in the fluid, on the average biomass. In particular, we considered a 2D kinematic velocity field that provides an idealized representation of the flow rolls observed in thermal convection [27], or in Langmuir circulation [START_REF] Craik | [END_REF][23][24] in the ocean mixed layer. Such flow pattern is characterized by both upwelling and downwelling regions, which play a key role on light-limited phytoplankton growth dynamics [4][5][6][7]. Indeed, a general feature emerging from previous studies (taking into account the vertical decrease of light intensity but not its horizontal variability) is that vertical mixing due to convective motions reduces the growth of phytoplankton, and may even lead to extinction events [5,6]. This is largely due to the accumulation of the population in the downwelling regions associated with persistent large-scale vortical structures. Consequently, phytoplankton is efficiently transported to depth, where reproduction is inhibited by the lack of light [5,6]. Irregular fluid motions on smaller time and length scales only weakly affect this phenomenology [6].

Our goal in this work was to explore if and how the response of the system, in terms of average biomass and survival conditions, changes when the incident light is not uniform at the surface, a situation that is often encountered in natural environments. A remarkable instance of such a configuration is found in oceanic waters that are partially covered by ice. In spite of their relevance for primary-productivity estimates [9,11], under-ice blooms are still not well understood. Notwithstanding its idealized character, the model considered here can then be informative about the basic mechanisms controlling the interplay between transport and biological growth processes in polar, light-limited environments, in the prespective of properly representing fluid motions in more complex numerical models.

Aiming to explore the simplest possible settings, we considered obstacles that fully absorb light and entirely cover either the downwelling region (LBD case) or the upwelling one (LBU case). We found that the main role of advective transport is to hinder phytoplankton growth. In both the LBD and LBU cases, the per-capita growth rate decreases with increasing flow intensity, and for sufficiently large advecting velocity, it can be become negative, meaning that the population goes extinct. While this behavior is qualitatively the same as the one observed in the absence of horizontal light-intensity modulation, now the quantitative effect depends on the position of the obstacles with respect to the regions of ascending and descending fluid motion. Phytoplankton survival is more difficult when the transmission of light is blocked over the upwelling flow region (LBU).

The flow field alone cannot be responsible for the asymmetric response of the system in the LBD and LBU cases, as it can be understood by examining the advection-reaction dynamics in a Lagrangian reference frame. We then explored the role of the sinking speed (due to the mass-density difference between phytoplankton organisms and water). Our analysis shows that, while typically much smaller than the flow intensity, the sinking speed is at the origin of the differences observed in the per-capita growth rates. We then provided a qualitative explanation of this phenomenon, which essentially depends on how the sinking speed combines with the vertical advecting velocity. In downwelling regions they sum, while in upwelling regions the sinking speed is subtracted from the flow. As a consequence, the residence time in the well-lit upper layer is increased in the upwelling region, relevant for the LBD case, and decreased in the downwelling region, relevant for the LBU case.

Finally, we have shown that the combination of velocities in the vertical discussed above is delicate, and can be strongly affected by lateral diffusion. Indeed, in the present system, vertical and horizontal motions cannot be easily decoupled. Diffusion spreads the population on the horizontal, bringing it from regions that are favorable for reproduction to deadly ones. This effect is more marked in the LBD case, where the population is more concentrated at the surface (in the upwelling region, where sinking partially counteracts the flow). In the LBU case, instead, the dynamics are less sensitive to the intensity of diffusion, because in the (downwelling) illuminated region the population is already transported (by both the flow and sinking) to deeper waters where reproduction is also strongly reduced.

To conclude, this study provided evidence that the onset of convection may be harmful to under-ice blooms, confirming previous findings [7]. The most notable result, however, concerns the role of the phytoplankton sinking speed. While small, the latter has a measurable impact on the growth dynamics of the population and can even be critical for its survival, in partially ice covered waters undergoing convective motions. This fact may have ecological relevance, as different phytoplankton species have different densities and, hence, different settling velocity. The system dynamics revealed highly complex, also in such an idealized setting, due to the very subtle interplay bewteen the different transport mechanisms (fluid advection, sinking, diffusion) and biological growth, which limits the possibility to obtain an analytical understanding. An interesting perspective would be to carry out numerical simulations in a more realistic configuration for the fluid, by considering the dynamics of turbulent convection in the presence of a melting interface with the ice obstacle [35].

Fig. 1 .

 1 Fig. 1. Schematic of the model setup. The vertical decrease of light intensity from the surface value Iin is described by Eq. (4).Phytoplankton (green dots) dynamics are assumed to evolve according to Eqs.[START_REF] Mann | Dynamics of marine ecosystems: biological-physical interactions in the oceans[END_REF] [with (3) and (4)] and (2), in the water columns in between the surface obstacles and below them (where Iin = 0), respectively. A cellular flow characterized by both upwelling and downwelling regions is also present.

Fig. 2 .

 2 Fig. 2. Instantaneous normalized population density field θ(x, z, t * )/ θ at a fixed instant of time t * = 600 h, for D = 10 cm 2 s -1 and U = (0, 0.71, 1.19) m h -1 [panels (a,d), (b,e) and (c,f), respectively], where θ stands for the spatial average. The obstacle is represented by the black rectangle, the vertical dashed lines delimit the light obstruction zone. The white lines are the isolines θ/ θ = 1, and t * is in the regime of stationary per-capita growth rate (rp(t) = rp = const). The solid black lines represent flow streamlines, with arrows indicating the circulation direction.

Fig. 3 .

 3 Fig. 3. Population per-capita growth rate rp(t), normalized by the intrinsic total (surface) growth rate r b , as a function of time, for various values of U , D = 10 cm 2 s -1 and Lz = 30 m. Panels (a) and (b) respectively refer to the LBD and LBU cases.

Fig. 4 .

 4 Fig. 4. Population per-capita growth rate, normalized by the intrinsic total (birth minus death) growth rate at the surface, as a function of time for different values of the sinking speed in both the LBD case and the LBU one. The flow intensity is U = 0.936 m/h, the diffusion coefficient is either D = 1 cm 2 s -1 (a) or D = 10 cm 2 s -1 (b).

Fig. 5 .

 5 Fig. 5. Distribution of N = 2000 particles at a time t > 0 for D = 1 cm 2 s -1 (red circles) and D = 10 cm 2 s -1 (blue squares) in the LBD (a) and LBU (b) cases. In each case, at time t = 0 all particles occupied a small area close to the top of the region receiving light (green circles).

Table 1 .

 1 Parameters of the biological dynamics.

	Parameter	Value	Meaning
	κ bg	0.2 m -1	Background turbidity
	κ	1.5 • 10

-11

m 2 cell -1 Specific light attenuation of phytoplankton pmax 0.04 h -1 Maximal specific production rate l 0.01 h -1 Specific loss rate H 30 µmol photons m -2 s -1 Half-saturation constant of light-limited growth Iin 350 µmol photons m -2 s -1 Incident light intensity v sink 0.04 m h -1

Acknowledgements

We are grateful to M. Cencini for interesting discussions about the role of diffusion.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.