domain: sde spi

text:
Sergio Chibbaro

Jair Reyes

Maurice Rossi

Alfredo Soldati

Francesco Zonta
email: francesco.zonta@tuwien.ac.at

Coherent structure modification by a shear acting at the surface of a turbulent open channel

Recent experiments in the river Seine have revealed the presence of persistent large rolls. This finding has motivated the present numerical study of a turbulent open channel flow. Our pseudo-spectral direct numerical simulations are aimed at understanding how the shear at the surface, for instance produced by an external wind, can change vortical structures in the bulk. Simulations are run at a Reynolds number in the range [5 × 10 3 ÷ 10 4 ], in line with similar laboratory experiments. We investigate how flow structures located near the bottom wall, near the surface or inside the bulk are modified by a surface shear. Statistical signatures are extracted through Fourier analysis of the simulated fields, and typical wavelengths of vortices or streaks are computed. Moreover, instantaneous fields are used to show the existence of upwelling and downwelling motions. These motions can be hindered by the presence of a large enough surface shear. This condition turns out to be necessary for the existence of streaks at the surface as well. Our investigation clarifies the dynamics of vortices in an open channel flow at moderate Reynolds number, indicating that unsteady vortex structures are indeed present, but the existence of long coherent rolls cannot be accounted for by the sole presence of shear at the surface. Numerical studies with a much longer domain and possibly at higher Reynolds numbers are needed to provide a firm answer to the question.

Introduction

The presence of coherent structures in turbulent flows has motivated a huge research effort in the last decades [START_REF] Robinson | Coherent motions in the turbulent boundary layer[END_REF][START_REF] Holmes | Turbulence, Coherent Structures (Dynamical Systems and Symmetry[END_REF][START_REF] Panton | Composite asymptotic expansions and scaling wall turbulence[END_REF]. While coherent structures are also present in isotropic flows [START_REF] She | Intermittant vortex structures in homogeneous isotropic turbulence[END_REF], we focus here on the properties of coherent vortices in the presence of physical boundaries, and in particular for the archetypal case of an open channel flows. Vortical structures can manifest either as a permanent secondary flow in a cross-stream direction or as an instantaneous pattern of the flow dynamics. In both cases, they are particularly important in environmental situations, since they may play a major role in sediment transport and dispersion [START_REF] Guazzelli | A physical introduction to suspension dynamics[END_REF].

Triggered by the influential work of Barenblatt [START_REF] Barenblatt | Scaling laws for fully developed turbulent shear flows. part 1. basic hypotheses and analysis[END_REF], a major effort has been made to understand better wall flows at moderate and high Re numbers. In such flows, experiments and numerical simulations revealed the existence of different near-wall vortical structures, spatially correlated and with different temporal persistency [START_REF] Jiménez | The structure of intense vorticity in isotropic turbulence[END_REF][START_REF] Adrian | Hairpin vortex organization in wall turbulence[END_REF]. Particularly interesting for boundary layers are the long streamwise elongated vortices [START_REF] Kline | The structure of turbulent boundary layers[END_REF][START_REF] Rogallo | Numerical simulation of turbulent flows[END_REF][START_REF] Jiménez | The autonomous cycle of near-wall turbulence[END_REF]. These streamwise structures seem to be 100-1000 long in terms of the wall unit distance [START_REF] Jiménez | The minimal flow unit in near-wall turbulence[END_REF], defined as the viscous scale δ ν ν/u τ where ν is the kinematic viscosity and friction velocity is u τ √ τ w /ρ, based on the averaged wall stress τ w and the fluid density ρ. Along the spanwise direction, these vortices generate streaks that generally extend over distances of approximately 100 wall unit distance or less. More recently, visualizations, experiments and numerical studies have revealed an even more complex situation with the existence of large-scale (LS) and very-large-scale motions (VLS) [START_REF] Kim | Very large-scale motion in the outer layer[END_REF][START_REF] Baltzer | Structural organization of large and very large scales in turbulent pipe flow simulation[END_REF][START_REF] Smits | High-Reynolds number wall turbulence[END_REF]. The typical length scale of these structures is the outer scale h rather than the viscous scale, and they can extend over very long distances in the streamwise direction (even up to 10h). It is worth noting that all these structures may be energetically important, but they are generated randomly in time and space, and are instantaneous representation of the Reynolds stress dynamics [START_REF] Pope | Turbulent flows[END_REF].

In hydraulic flows, such as open channels, the situation is less clear-cut. In rough open channels, there is some evidence that large structures are present as streamwise rotating rolls, which actively participate to the bed dynamics [START_REF] Barros | Observations of turbulent secondary flows in a rough-wall boundary layer[END_REF][START_REF] Franca | Detection and reconstruction of large-scale coherent flow structures in gravel-bed rivers[END_REF][START_REF] Cameron | Very-large-scale motions in rough-bed open-channel flow[END_REF][START_REF] Vowinckel | Spatially-averaged momentum fluxes and stresses in flows over mobile granular beds: a dns-based study[END_REF]. Smooth channels are equally relevant when the interaction with the rough bed may be neglected. In this case, the presence of streaky structures at the bottom has been corroborated by many studies [START_REF] Gulliver | Measurements of large streamwise vortices in an open-channel flow[END_REF][START_REF] Shvidchenko | Critical shear stress for incipient motion of sand/gravel streambeds[END_REF]. Similar low-and high-speed streaks have also been found in free surface flows [START_REF] Rashidi | Turbulence structure in free-surface channel flows[END_REF][START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF][START_REF] Kumar | An experimental investigation of the characteristics of free-surface turbulence in channel flow[END_REF]. Streaks exist at the free surface at least when some large enough external shear is imposed at the free surface. This result is non-trivial since, contrary to a solid wall, vortex lines can attach at a free surface and a wall-normal vorticity component can exist. This is one classical instance of modification of vorticity dynamics by a free surface [START_REF] Sarpkaya | Vorticity, free surface, and surfactants[END_REF]. Similarly to what was proposed for closed wall flows, some authors have linked the presence of streaks near the surface to vortical motions spanning the entire channel height in the cross-stream direction rather than being localized near boundaries [START_REF] Kaftori | Funnel-shaped vortical structures in wall turbulence[END_REF][START_REF] Tamburrino | Large flow structures in a turbulent open channel flow[END_REF][START_REF] Tamburrino | Free-surface visualization of streamwise vortices in a channel flow[END_REF][START_REF] Adrian | Coherent structures in flow over hydraulic engineering surfaces[END_REF][START_REF] Sukhodolov | Flow dynamics in alluvial channels: the legacy of Kirill v[END_REF][START_REF] Zhong | Coherent structures and their interactions in smooth open channel flows[END_REF].

In addition, in smooth open channel flows, some in situ experiments have pointed out the presence of vortices, which look like secondary flows [START_REF] Nikora | Secondary flows in rivers: theoretical framework, recent advances, and current challenges. Gravel-bed rivers: Processes, tools, environments[END_REF] in the transverse plane, i.e., contrary to LS and VLS in pipe and boundary layer flows, they are not random in space and time. Under which conditions these vortices are encountered in natural flows, and what is their impact on sediment/species transport, hydraulic resistance, mixing and morphodynamics are still a matter of debate that requires further analyses [START_REF] Cameron | Very-large-scale motions in rough-bed open-channel flow[END_REF][START_REF] Nikora | Secondary flows in rivers: theoretical framework, recent advances, and current challenges. Gravel-bed rivers: Processes, tools, environments[END_REF].

To provide some insights on this issue, in the present work, we study numerically an open channel flow, focusing in particular on the upwelling and downwelling motions of fluid, and investigating whether there exists a clear statistical signature of rolls scaling with the height of the channel. Some recent works have particularly motivated the present one. Experiments by Chauvet et al. [START_REF] Chauvet | Recirculation cells in a wide channel[END_REF] have analyzed velocity and vorticity in the river Seine. Through Doppler measurements, these authors have provided evidence of a secondary flow in the transverse plane. More precisely they have measured a time-averaged vertical velocity, thus showing the existence of a typical wavelength along the spanwise direction that may be related to the presence of large circulation rolls. These structures appear to be robust, yet they are of much smaller amplitude than the main streamwise component: For the river Seine, the average speed is of about 1.6 [ms -1 ], whereas the vertical velocity turns out to be of the order of 0.510 -2 ms -1 . How these structures are formed is not yet clear. Of even more interest for the present work, Zhong et al. [START_REF] Zhong | Statistical analysis of turbulent super-streamwise vortices based on observations of streaky structures near the free surface in the smooth open channel flow[END_REF] recently studied the dynamics of large vortical structures performing laboratory experiments in a smooth and large open channel (i.e., with a large width-to-height aspect), at Reynolds numbers similar to those used in our study. Given the close similarity with our study, we will make explicit reference to Zhong et al. [START_REF] Zhong | Statistical analysis of turbulent super-streamwise vortices based on observations of streaky structures near the free surface in the smooth open channel flow[END_REF] in discussing our results.

In this context, as pioneering shown by Lam and Banerjee [START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF] and Pan and Banerjee [START_REF] Pan | Numerical simulation of particle interactions with wall turbulence[END_REF], direct numerical simulations (DNS)-because of the large amount of data and information that generate-can be a valuable tool to probe the physics of turbulent open channel flow and to test hypotheses and theories. The major disadvantage of DNS is the limited value of the Reynolds number that can be achieved, which is typically rather small compared to that of geophysical flows. This has motivated the use of large-eddy simulations (LES) [START_REF] Calmet | Statistical structure of high-reynolds-number turbulence close to the free surface of an open-channel flow[END_REF][START_REF] Rodi | Large-eddy simulation in hydraulics[END_REF] for large-scale geophysical flow simulations, despite the tendency of LES to smoothen out streaky structures [START_REF] Bianco | Intrinsic filtering errors of Lagrangian particle tracking in les flow fields[END_REF].

The main goal of the present work is thus to provide a careful statistical analysis of a turbulent open channel flow, in order to understand whether large vortices appear and whether their manifestation may be linked to the presence of streaks at the surface. For such a purpose, we decided to run DNS. The setup and the value of the main parameters of the simulations are consistent with those of similar laboratory experiments [START_REF] Zhong | Coherent structures and their interactions in smooth open channel flows[END_REF][START_REF] Zhong | Statistical analysis of turbulent super-streamwise vortices based on observations of streaky structures near the free surface in the smooth open channel flow[END_REF]. Because of the moderate Reynolds number considered here, the work is intended firstly as a complement to these ones. Furthermore, as in 2D simulations at low Re [START_REF] Borue | Interaction of surface waves with turbulence: direct numerical simulations of turbulent open-channel flow[END_REF], several shear rates are imposed at the surface via the boundary conditions to assess the possible effect of the shear stress in triggering the growth of large vortices. As mentioned above, the ratio between the amplitude of the rolls and the mean flow is quite small in natural conditions, and therefore, we have performed simulations at high external shear compared to natural configurations, to improve the signal/noise ratio in the hope of better revealing such structures.

Methods

Let us consider a turbulent open channel flow driven by a mean pressure gradient, which is formerly equivalent to study the case of a free surface flow flowing by gravity over an inclined plane. The geometric configuration is represented by a rectangular domain of length Lx along the streamwise direction, height Lz and width L y . We assume: Lx 4π ĥ, L y 2π ĥ and Lz 2 ĥ,

where ĥ is the half channel height (Fig. 1), symbol • indicating dimensional quantities. A Cartesian frame is used, where the streamwise, spanwise and normal coordinates are denoted, respectively, x, ŷ and ẑ. The flow is assumed to be incompressible:

∂ û j ∂ x j 0, (2)
and governed by the Navier-Stokes equations

∂ ûi ∂ t + û j ∂ ûi ∂ x j - 1 ρ ∂ p ∂ xi + ν ∂ 2 ûi ∂ x j ∂ x j , ( 3 )
where û ≡ û1 , v ≡ û2 and ŵ ≡ û3 are the streamwise, spanwise and normal velocity components and p pressure. Fluid density ρ and kinematic viscosity ν are assumed to be constant. Periodic conditions are imposed along the spanwise direction ŷ, so to mimic a flow that is contained in a channel but located away from lateral walls and therefore not much influenced by them. Along the streamwise direction x, periodicity is also imposed for the velocity field. For pressure, we assume that there exists a constant mean pressure gradient ˆ and that pressure perturbations around the mean value are periodic along the streamwise x and spanwise L y and height Lz . The origin is located at ( x, ŷ, ẑ) (0, 0, 0) so that bottom and surface are situated, respectively, at ẑ -ĥ and ẑ ĥ, and horizontal coordinates x, ŷ are positive ŷ directions. The mean pressure is related here to boundary conditions, but it could be also introduced from constant forces in Eq. ( 3) representing the gravity component along the bed. A larger length is imposed along x to take into account the correlation of streamwise vortices [START_REF] Jiménez | The minimal flow unit in near-wall turbulence[END_REF]. In our simulations, a no-slip condition is imposed at the bottom, i.e., ( û, v, ŵ) (0, 0, 0) at ( x, ŷ, -ĥ) and the free surface displacements along the vertical direction are assumed to be small with respect to the channel depth. This implies that the vertical velocity component and the vertical surface displacements can be neglected: ŵ 0 at ẑ ĥ. This rigid-lid approximation eliminates surface waves. While the presence of waves may add important effects [START_REF] Morland | Waves generated by shear layer instabilities[END_REF][START_REF] Ambrosi | Non-normal stability analysis of a shear current under surface gravity waves[END_REF], this choice is considered as an adequate first order model for the purpose of the present work [START_REF] Lovecchio | Influence of thermal stratification on the surfacing and clustering of floaters in free surface turbulence[END_REF][START_REF] Zonta | Growth and spectra of gravity-capillary waves in countercurrent air/water turbulent flow[END_REF], which considers flow at moderate Reynolds numbers.

For streamwise and spanwise velocity components, we impose at the surface ẑ ĥ either a free-slip condition (FSC) or an imposed stress condition (ISC). The FSC corresponds to a free surface with no imposed shear (that is, neglecting the shear imposed by the atmosphere). This leads to the following conditions at the free surface:

μ ∂ ∂ ẑ û( x, ŷ, + ĥ) 0 and μ ∂ ∂ ẑ v( x, ŷ, + ĥ) 0, (4)
where μ ρ ν is the dynamic viscosity. For the ISC, we impose a shear τs 0 at the free surface plane ẑ ĥ which is directed along the stream and results in the following conditions at the free surface:

μ ∂ ∂ ẑ û( x, ŷ, + ĥ) τs and μ ∂ ∂ ẑ v( x, ŷ, + ĥ) 0. (5)
Because of the presence of a mean pressure gradient and of an imposed shear at the surface, this flow may be called a turbulent Poiseuille-Couette flow. For the FSC, the global equilibrium of forces between mean pressure gradient ˆ and the mean shear stress at the bottom wall τw gives τw -2 ĥ ˆ . For the ISC, the mean pressure gradient is always balanced by a combination of the shear at the bottom τw and the shear at the surface τs , that is, τw -τs -2 ĥ ˆ .

The system is solved in dimensionless form by introducing a characteristic length [L] and velocity [U]. The length [L] is based on the half channel height ĥ and the velocity on the mean pressure gradient, i.e., [U ] ĥ ρ | ˆ |. The following dimensionless variables are thus defined: (x, y, z)

1 [L] ( x, ŷ, ẑ), (u, v, w) 1 [U ] ( û, v, ŵ), t 1 [L][U ] -1 t, p 1 ρ[U ] 2 p, τ 1 ρ[U ] 2 τ .
In dimensionless form, the computational domain becomes x : [0, 4π], y : [0, 2π] and z : [-1, +1], and the governing Eqs. ( 2) and (3) read:

∂u j ∂ x j 0, ( 6
)
∂u i ∂t + u j ∂u i ∂ x j - ∂ p ∂ x i + 1 Re ∂ 2 u i ∂ x j ∂ x j . ( 7
)
where the Reynolds number is defined Re ≡ ĥ ĥ From a numerical point of view, Eqs. ( 6) and ( 7) are solved via a standard pseudo-spectral method [START_REF] Canuto | Spectral methods[END_REF][START_REF] Zonta | Modulation of turbulence in forced convection by temperature-dependent viscosity[END_REF]. Further details on the numerics can be found in Appendix A. Table 2 Flow conditions in dimensional units: ĥ is the height, ν the kinematic viscosity, Ŵ the channel width, Ŵ / ĥ the aspect ratio. To make the comparison, we have used the definitions from Zhong et al. [START_REF] Zhong | Statistical analysis of turbulent super-streamwise vortices based on observations of streaky structures near the free surface in the smooth open channel flow[END_REF].Note that, for this reason, the Re τ here corresponds to the double of that obtained in We present five numerical simulations of the turbulent Poiseuille-Couette flow: a single FSC and four ISC simulations. In each ISC simulation, the imposed shear stress τ s corresponds to a given mean streamwise velocity V s at the surface. These simulations are sorted in ascending order according to the imposed shear stress τ s . Simulation ISC 0 with the imposed velocity value V s 0 is similar, but not identical, to a standard Poiseuille flow since, differently for the standard Poiseuille flow, the stream-and spanwise velocity components can fluctuate at the top surface. ISC 0 and ISC 1 correspond to a negative shear τ s while ISC 2 and ISC 3 correspond to a positive shear τ s . The different parameters are presented in Table 1. In order to indicate which real experiment could be directly related to our results, we present in Table 2 the case without shear applied at the surface with typical dimensional units, together with the parameters of the experiments reported in [START_REF] Zhong | Statistical analysis of turbulent super-streamwise vortices based on observations of streaky structures near the free surface in the smooth open channel flow[END_REF].

To run the above simulations, we proceed as follows. We start from a realization of a closed channel flow, we changed the boundary condition at the top wall from a no-slip boundary into a free-slip one (FSC), and we let this new simulation to evolve until a statistically steady state is reached. We then use a single realization of this field as initial condition for the other simulations. For each simulated case, after an initial transient, the flow attained a new statistically steady state, in which statistics can be computed. All the statistics that will be shown below have been checked to be converged with respect to the number of mesh points and time discretization.

Results

Analysis of velocity statistics

Since we are particularly interested in coherent structures, we decided to present most of our results concerning classical statistics in appendix (see appendix B-C). As a brief statement, current results confirm previous analyses carried out at lower Reynolds numbers [START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF][START_REF] Pan | Numerical simulation of particle interactions with wall turbulence[END_REF]. For instance, we clearly observe (see Table 1) that, compared to FSC, the mean streamwise velocity u m decreases for ISC 0 and ISC 1 , and increases for ISC 2 and ISC 3 . The boundary condition at the top surface influences the local gradient of u near the bottom. Compared to FSC, we find larger (respectively, smaller) bottom shear stress when the imposed mean surface velocity V s is larger (respectively, smaller) than the mean surface velocity of FSC. The wall-normal profile of the dimensionless mean shear stress, given by:

τ (z) 1 Re d u dz -u w , ( 9 )
with velocity fluctuations u i u iu i , is expected to be τ (z) τ s + (1z) and τ w τ s + 2, where τ w ≡ τ (-1) and τ s ≡ τ (1). The agreement of the numerical shear stress with the profile given by Eq. ( 9) is excellent (see appendix B). This result constitutes also a check of the statistical convergence of our computations. To properly quantify the effects due to viscosity near the two boundaries, we can further analyze the mean stress by decomposing it into its viscous and turbulent part and by rescaling the viscous and turbulent stresses as follows. For the bottom boundary, the distance from the bottom is measured in wall units as

z + ≡ 1 δ ν (z + 1) with δ ν ≡ 1 Re τ , ( 1 0 )
with δ ν the characteristic length scale of the viscous sublayer. Quantity z + can be viewed as a local Reynolds number based on velocity u τ and distance from the wall (z + 1). In the same manner, we measure the distance from the surface in surface units as

z + s ≡ 1 δ ν s (1 -z), with δ ν s ≡ 1 Re τ s ( 11
)
with δ ν s the viscous sublayer length scale at the surface. In Fig. 2, the ratio Re -1 τ (z) d u dz (respectively, -1 τ (z) u w ) of the viscous stress (respectively, Reynolds stress ) to mean shear stress τ (z) is plotted with emphasis on the viscous layer near the bottom (as a function of z + ) and near the surface (as a function of z + s ). Near the bottom boundary, these profiles nicely collapse regardless of the imposed τ s . A similar behavior is found near the surface boundary for the different ISC, with the possible exception for ISC 2 , which corresponds to the smallest imposed shear τ s . ISC 2 is probably a limit case where the surface does not differ much from a solid wall. For FSC, the structure of shear is different near the surface since both viscous and Reynolds stress are zero at the surface, (Hence, it is not shown in Fig. 2b.)

These results together with the statistical analysis of mean streamwise velocity and fluctuations (see Appendix B-C) allow us to conclude the following: (i) Considering previous results obtained at lower Re number, no appreciable Reynolds number effect is found on the main statistics of the flow, at least for the range of Reynolds number considered here; (ii) the flow near the surface is similar to a boundary layer, as far as the mean velocity is concerned; (iii) the same conclusion applies to fluctuations statistics; and (iv) the stresses, upon proper rescaling, appear independent of the imposed external shear.

Flow structures at the bottom boundary

We switch now to the main focus of this work, namely coherent flow structures. Due to the significant shear stress present at the bottom boundary, the near-wall region is dominated by local streamwise vortices and low-speed streaks (velocity deficits). These structures have also been widely observed in laboratory experiments [START_REF] Rashidi | The effect of boundary conditions and shear rate on streak formation and breakdown in turbulent channel flows[END_REF][START_REF] Babiker | Vortex imprints on a free surface as proxy for surface divergence[END_REF] and reproduced in numerical simulations at low Reynolds number [START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF]. To characterize the near-wall dynamics in the present problem, in Fig. 3 we plot the contour maps of the fluctuations u of the streamwise velocity on the plane z + 5 at the upper limit of the viscous sublayer. Results are shown for the cases ISC 0 , FSC and ISC 3 only, and only a part of the domain is shown to highlight the structures (in appendix D, the plots of all cases and in the entire domain can be found). Figure 3, as well as the subsequent plots displaying instantaneous fields, presents the general behavior. The flow has similar characteristics for all cases, with the presence of low-and high-speed streaks populating the near-wall region. However, the length and width of these streaks do depend on the applied shear. In particular, the width of streaks decreases with decreasing the bottom shear. A more quantitative indication of the streak separation is obtained from the average Fourier spectrum of the vertical velocity component w as a function of the spanwise wavelength λ y . This quantity is computed at z + 5 and averaged in space (along the x direction) and in time (Fig. 4).

For all cases, the spectrum displays a peak that represents a characteristic wavelength λ y : 0.53 for ISC 0 , 0.45 for ISC 1 , 0.39 for FSC, 0.33 for ISC 2 , 0.28 for ISC 3 . Streaks exhibit a shorter separation as the shear stress at the bottom boundary is increased (which is in turn a consequence of the increased module of the shear at the top surface). Nevertheless the same quantity normalized by the shear velocity w /u τ and presented against the wavelength expressed in wall units λ + y ≡ 1 δ ν λ y , displays similar characteristic separation λ + y (z + 5) ≈ 90 ÷ 100 for all cases, see also Fig. 6. To get an idea on how these patterns are changing with respect to the distance from the wall, two pictures are presented. First (Fig. 5), the instantaneous streamwise velocity u is plotted in the transverse plane yz. Low-velocity streaks, having a height of about 40-50 wall units, are observed for all cases. Second (Fig. 6), the average Fourier spectrum of w /u τ is also presented at three positions z + within the buffer layer (z + < 30), as a function of wavelength expressed in wall units λ + of streaks does not depend much on the different cases at any particular position z + . However, streaks spacing slightly increases with the distance from the wall: λ + y (z + 5) ≈ 90, λ + y (z + 15) 100, λ + y (z + 30) 135. The λ y values for the FSC case agree with those obtained by Kim et al. [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF], in particular λ + y ≈ 100 for z + 10. Beyond the buffer layer, viscosity effects are negligible and wall units cannot adequately characterize the flow. Hence, the dominant wavelengths λ + y computed from the different cases do not collapse when considering planes at z + 50. By contrast, for such planes, the wavelength at which the maximum amplitude is reached is similar in terms of λ y for the different cases (Fig. 7). Furthermore, for z + > 50 spectra are similar and amplitudes do not change much. This confirms that streaks have a height of the order of 50 wall units.

Flow structures at the free surface

One issue concerns the interaction of turbulence with the free surface and the effect of the shear imposed at the surface on the flow structures close to the surface. In Fig. 8 (and in appendix D), we present the fluctuations u of the streamwise velocity component in the horizontal plane located at z + s 5 (in surface units) (but for FSC, for which the plot is at 1-z δ ν 5). Streak-like patterns are clearly visible for ISC 0 and ISC 3 (see Fig. 8 top and bottom, respectively), and also ISC 1 . For FSC (see Fig. 8 center) and ISC 2 , a streak-like pattern is not observed. This suggests that flow streaks are likely to appear only when the shear exceeds a critical value, as previously found by [START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF]. To quantify these observations, we compute (i) the average Fourier spectrum of the vertical velocity component as a function of the spanwise wavelength λ y (see Fig. wavelength corresponds to a value λ + y, s 80, similar to the value obtained for the bottom streaks. As a consequence, elongated structures of streamwise vorticity appear at the surface for ISC 0 , ISC 1 and ISC 3 , but not for FSC and ISC 2 .

Horizontal velocity divergence at the surface

To characterize further the interaction between the flow structures in the bulk of the channel and the free surface, we compute the horizontal flow divergence:

∇ H • u ≡ ∂u ∂ x + ∂v ∂ y ( 12
)
at the surface, i.e., at z 1. The two-dimensional flow divergence is associated with the exchange of mass between the surface and the bulk of the flow [START_REF] Lovecchio | Upscale energy transfer and flow topology in free-surface turbulence[END_REF]. Regions of positive horizontal divergence are regions of local flow expansion. These are due to incoming 10. When the shear at the surface is weak (ISC 2 ) or vanishing (FSC), upwellings (red regions) and downwellings (blue regions) are lumps of positive or negative ∇ H • u with no specific spatial orientation. For ISC 0 , ISC 1 (not shown for brevity) and ISC 3 in which a large shear is present, the situation changes drastically: Upwellings and downwellings appear as longitudinal structures. These motions are hence distorted by the shear applied at the surface and tend to align with high-and low-speed streaks. They are also weaker than for ISC 2 (not shown for brevity) or FSC. Therefore, the presence of a shear is found to disrupt upwellings and downwellings (large-scale flow events). This in turn allows the presence of streaks at the surface. To fully characterize the flow structure near the upper surface, the concept of the flow compressibility is introduced. A dimensionless compressibility factor [START_REF] Gawedzki | Phase transition in the passive scalar advection[END_REF], linked to the local flow divergence, can be computed on xy planes at different z locations:

C ( ∇ H • u ) 2 ∂u ∂ x 2 + ∂u ∂ y 2 + ∂v ∂ x 2 + ∂v ∂ y 2 . ( 13
)
This factor is nonzero because of the 3D incompressibility of the flow. Figure 11 shows its profile for the top half of the domain.

The free-shear case reproduces the result recently discussed by [START_REF] Lovecchio | Upscale energy transfer and flow topology in free-surface turbulence[END_REF]. The value observed at the surface (z 1) in these simulations is C 0.6, which is higher than that found in experiments and numerical simulations considering the free-shear surface of an otherwise 3D incompressible turbulent isotropic flow [START_REF] Cressman | Eulerian and Lagrangian studies in surface flow turbulence[END_REF]. In this case, the compressibility factor turns out to be C ≈ 0.5 [START_REF] Cressman | Eulerian and Lagrangian studies in surface flow turbulence[END_REF]. The value at the surface is close but higher also than the theoretical value C 0.5 that characterizes the strong compressible Kraichnan flow [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF][START_REF] Kraichnan | Anomalous scaling of a randomly advected passive scalar[END_REF]. This indicates an important compressibility effect in the surface case, even though it seems that time correlations of real flows tend to weaken the compressibility effect with respect to the Kraichnan model [START_REF] Boffetta | Lagrangian tracers on a surface flow: the role of time correlations[END_REF]. However, this is still a matter of debate, and the role of flow time correlations on the flow compressibility, and on the corresponding flow mixing, should be clarified better [START_REF] Dhanagare | Weak-strong clustering transition in renewing compressible flows[END_REF].

The compressibility factor decreases to a value near C ≈ 0.2 in the bulk of the channel. This value is close to the value 1/6 which is expected for a two-dimensional cut of a three-dimensional homogeneous, isotropic flow. The imposition of a shear (both positive or negative) does not modify the value of C in the channel core (z 0) but strongly reduces its value at the surface. Moreover, the value is decreased also in a large region below the surface 0.5 < z < 0.9. This should be related to the distortion of the upwelling and downwelling motions operated by shear. The ISC 2 case is the only one in which the effect is smaller. In this case, the value of C is almost everywhere 0.2 as in the bulk. The results show a two dimensionalization of the flow in the region below the surface in the presence of shear. Since streaks are quasi-2D structures, this finding is consistent with the fact that streaks are only observed when a sufficiently high shear is present (otherwise they are not observed).

Flow structures inside the channel

We investigate now whether a shear applied at the surface promotes or hinders the formation of large-scale structures inside the bulk. Away from the boundaries, the formation of recirculation cells may be identified by using the vertical velocity w in the transversal 13 Average Fourier spectrum of the vertical velocity component normalized by friction velocity as a function of the spanwise wavelength λ y at z -0.5 (left panel), the center of the channel, z 0, (central panel) and at z 0.5 (right panel) yz plane (Fig. 12, left column) and in horizontal xy plane (Fig. 12, right column). Some difference between cases can be pointed out. For FSC (and ISC 2 , not shown here), large-scale structures generated in the central region of the channel can extend up to the surface (central panel in Fig. 12). They consist of adjacent upwelling and downwelling structures alternating in space. Without a shear (or with a small positive one), these large-scale vertical motions occur with a rather random horizontal spatial distribution, as highlighted in Fig. 12 (right column). A negative shear at the surface of the channel promotes the formation of smaller rollers close to the surface, with the consequent weakening of the large-scale structures (see case ISC 0 , top panel in Fig. 12). The results obtained in the ISC 3 case (bottom panel in Fig. 12) might suggest that a strong positive shear increases the persistence along the stream of large-scale rolls.

To quantify these observations, we compute the average Fourier spectrum (on a xy plane) of the vertical velocity component as a function of the spanwise wavelength λ y . Results are shown in Fig. 13 for three different vertical positions (z -0.5, z 0 and z 0.5). Regardless of the vertical location of planes, the dominant wavelength for the case of surface positive shear (ISC 2 and ISC 3 ) is at λ y ∼ 2, i.e., the channel height. In addition, a positive shear seems to increase the amplitude of the large-scale rolls. For negative shear cases (ISC 0 and ISC 1 ), the spectrum shows that the amplitude is reduced compared to the previous two cases and the optimal length is also less pronounced as we approach the wall or the free surface (z -0.5, z 0.5). These evidences corroborate the existence of some large-scale coherent dynamics, at least of modest intensity because of the relative importance of the peak in the Fourier spectrum. In order to highlight whether the vertical movements are correlated in the horizontal direction, the vertical velocity is visualized in the horizontal plane at the center of the channel (Fig. 12 right column). From these results, it is possible to argue that large-scale structures might be also more correlated when a sufficient and positive shear is imposed at the surface (ISC 3 ), indicating that a strong favorable wind triggers the formation of large-scale vortices, as in similar hydrodynamic unstable shear flows, i.e., Couette flows. Yet, the streamwise correlation remains poor, and results cannot be considered conclusive.

Discussion and conclusions

We have carried out a numerical analysis of the dynamics of vortical structures in turbulent open channel flows. This issue, which is of great fundamental interest, is also relevant for river dynamics, for which experimental and numerical studies seem to highlight the presence of large elongated rolls, whose size in the cross-plane scales with the river height h, and whose length can be several times h. Our numerical experiments have been made with parameters that can be considered similar to the recent laboratory experiments by Zhong et al. [START_REF] Zhong | Statistical analysis of turbulent super-streamwise vortices based on observations of streaky structures near the free surface in the smooth open channel flow[END_REF]. This numerical study complements those laboratory experiments, since numerical simulations have the advantage of providing all possible information with great accuracy. We have carried out different simulations choosing different shears applied at the surface, which mimic the presence of a wind, either opposite or along the streamwise direction. It is worth emphasizing that our simulations are made in idealized conditions and comparison with field experiments, always at much higher Reynolds numbers, must be taken cum grano salis. We discuss here our main points.

(i) Streaks at the surface are found only when a sufficient shear is imposed at the surface, and the separation of these streaks is found to be of the order of 50 ÷ 100δ ν . In this context, our results confirm and extend previous results obtained at lower Re [START_REF] Pan | Numerical simulation of particle interactions with wall turbulence[END_REF], hence indicating, at least for the considered range of Reynolds numbers, very small (if any) Reynolds number effect. For the range of parameters investigated here, the observed streamwise vortices seem to scale in inner units, as those encountered in near-walls flows.

(ii) The boundary layer at the bottom is almost unaffected by the shear at the surface. Therefore, an external shear does not seem to be responsible for important changes in the near-wall region (river bed).

(iii) We have investigated the dynamics of upwelling and downwelling motions, and the influence of the shear on them. Upwelling and downwelling motions are always present along the full height of the channel, but the applied shear can modify their dynamics:

123 While in the free shear or in the presence of a very low one, these motions are randomly organized in the horizontal plane, in the other cases, these motions are more regular and oriented, and streaks are also more visible.

(iv) Our analysis indicates that the streaky structures, and in particular their size, are sensitive to varying the amplitude and the direction of the shear. We have also tried to understand whether the structures, present along the vertical direction, extend also for several channel height in the streamwise one. A strong positive shear seems to trigger the horizontal organization of large rolls. Yet this result cannot explain the persistency of such vortices found in several in situ experiments. When a negative shear is present, the vertical structures appear not to be well correlated in the streamwise direction. That seems to imply that a low negative shear is probably necessary to deform the upwelling and downwelling movements, but is not sufficient to justify the presence of large-scale vortices. It is worth noting that the general case for in situ experimental observations would be of negative shear. In any case, we have not found in our simulations sufficient signature of large elongated roll in smooth open channels due to shear. Those large elongated rolls found by other authors are not produced by simple shear at the surface, at least in the conditions considered here. Large rolls are in any case of small amplitude, an occurrence that makes their observation somehow tricky. Our results are in agreement with recent laboratory experiments [START_REF] Zhong | Coherent structures and their interactions in smooth open channel flows[END_REF][START_REF] Zhong | Statistical analysis of turbulent super-streamwise vortices based on observations of streaky structures near the free surface in the smooth open channel flow[END_REF], in which authors measured a small correlation in the vertical and longitudinal plane, and linked it to large streamwise vortices. The experiments were performed at similar Reynolds numbers, but the domain was longer, possibly allowing for the development of larger rolls. Our analysis of the correlation in the horizontal plane suggests that these results cannot be considered conclusive for the existence of coherent large rolls.

According to our findings, the presence of a shear seems a necessary, but not sufficient, condition. Hence, large-scale structures should arise via some other physical mechanisms. The first idea is to invoke the effect of very high Re number [START_REF] Smits | High-Reynolds number wall turbulence[END_REF]. Unfortunately, a parametric study of a well resolved long open channel at high Reynolds is prohibitive from the computational point of view. Stratification may trigger also large convection rolls, but in hydraulic flows it would be stable in most of the cases, making this hypothesis improbable. A promising perspective is to look at the possible effect of surface waves, which are known to be responsible for large vortices, as the Langmuir circulation in oceans [START_REF] Leibovich | The form and dynamics of Langmuir circulations[END_REF][START_REF] Thorpe | Langmuir circulation[END_REF]. We plan to make some numerical simulations with a moving free surface, instead of using the rigid-lid approximation to analyze the effect of waves on the fluid motion.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Numerical method

The Navier-Stokes equations are solved applying the no-slip condition at the bottom, (u, v, w) (0, 0, 0) at z -1. At the top surface z 1, the rigid-lid condition becomes w 0 and ∂ ∂z u(x, y, + 1) 0, ∂ ∂z v(x, y, + 1) 0, for FSC, ( 14)

∂ ∂z u(x, y, + 1) Re τ s , ∂ ∂z v(x, y, + 1) 0, for ISC. ( 15
)
This set of governing equations is discretized using a pseudo-spectral method [START_REF] Canuto | Spectral methods[END_REF][START_REF] Zonta | Modulation of turbulence in forced convection by temperature-dependent viscosity[END_REF], which is based on transforming the field variables into wavenumber space, through Fourier representations for the (homogeneous) periodic directions x and y and Chebyshev polynomials T n z (z) cos(n z arccos(z)) for the non-homogeneous wall-normal direction z u i (x, y, z; t) ≡

|n x | |n y | |n z | ũi (k x , k y , n z ; t)T n z (z) e i(k x x+k y y) . ( 16
)
The summation is performed over the integer n x , n y and n z , which vary in the range -N x /2 + 1 ≤ n x ≤ N x /2, -N y /2 + 1 ≤ n y ≤ N y /2 and 0 ≤ n z ≤ N z where N x , N y and N z -1 are powers of two. The associated spatial grid points are, respectively, N x , N y and N z along the x, y and z directions

x(i) (i -1) L x N x -1 i 1, 2, 3, ..., N x (17a)
y( j) ( j -1) L y N y -1 j 1, 2, 3, ..., N y (17b) z(k) cos k -1 N z -1 π k 1, 2, 3, ..., N z . ( 17c
)
Time advancement of the governing equations is achieved by a combination of an implicit Crank-Nicolson scheme for the viscous terms and an explicit Adams-Bashforth scheme for the convective nonlinear terms. As commonly done in pseudo-spectral methods, the convective nonlinear terms are first computed in physical space and then transformed in the wavenumber space using a de-aliasing procedure based on the 2/3 rule [START_REF] Canuto | Spectral methods[END_REF]; derivatives are evaluated directly in the wavenumber space to maintain spectral accuracy. Further details about the numerical approach can be found in the literature [START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF][START_REF] Zonta | Modulation of turbulence in forced convection by temperature-dependent viscosity[END_REF]. As far as computational time is concerned, each of the 5 simulations has taken about 4 weeks using 256 cores on an Intel PC cluster. A summary of the main simulation parameters is given in Table 3

B Mean velocity and stress

In Fig. 14a, we present the wall-normal behavior of the mean streamwise velocity profile u for all cases. Although the behavior of the mean streamwise velocity for the case FSC is similar to that observed in one-half of a turbulent Poiseuille channel flow, there exists a difference between the two cases: Vorticity is normal to the surface in FSC but not at the center of a turbulent Poiseuille flow.

For ISC, we observe (Fig. 14c) that the top boundary condition enables momentum transfer by diffusion at the free surface (boundary layer behavior close to the free surface). At the free surface, the rigid-lid condition imposes a zero Reynolds stress (Fig. 14d) through the prescription of w 0. The mean streamwise velocity profile (see Fig. 15) is computed in wall units as:

u + ≡ u u τ Re Re τ u . ( 18
)
For all the cases, within the viscous sublayer (z + < 5), u + collapses onto the linear behavior u + z + prescribed by the law of the wall. In the core region (z + > 30), the logarithmic law is satisfied up to 40% of the half channel height h, although with slightly different values for the von Kármán constant κ and translation constant C for the various cases, as verified through data fitting. The behavior is similar to what was found previously [START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF] at lower Reynolds number. In addition, there is no clear trend for the values of κ and C for increasing u τ . For FSC, we found κ 2.5 -1 and C 4.94: The value of κ is the same obtained by [START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF], and their constant C 5.1 is very slightly larger. That is probably within the combined error but it may also be due to the fact that present simulations are at a higher Re.

We also characterize the mean streamwise velocity profile close to the free surface when a shear stress is applied (ISC cases). To do this, the surface mean velocity u(1) ≡ F s is subtracted to the total mean velocity profile. The velocity given by:

u + s ≡ abs( u -V s ) u τ s . ( 19
)
which is measured in surface units (Fig. 16). The velocity profiles follow the law of the wall with characteristic constants that depend on the flow configuration. This means that the velocity profile can be approximated in the viscous sublayer, respectively, in the core by

u + s z + s respectively u + s 1 κ ln z + s + C. ( 20
)
The viscous sublayer at the surface (0 ≤ z + s 3) is smaller than its counterpart at the bottom. Similarly, the characteristic coefficients 1 κ and C are generally smaller than their counterparts at the bottom, but no clear trend is observed. From the above considerations, we conclude that the flow near the surface is similar to a boundary layer as far as the mean velocity is concerned.

C Statistics of velocity fluctuations

Regarding fluctuations statistics, we focus on the root mean square (RMS) u i rms ≡ u i 2 of the velocity component u i . The wall-normal behavior of u i rms for the three fluctuating velocity components normalized by the friction velocity u τ is shown in Fig. 17 for a region close to the bottom wall 0 < z + < 80 (the upper limit being about z ≈ 0.4). As commonly observed in wall-bounded turbulence, the values of u rms are larger than the corresponding values of v rms and w rms . For all cases, the u rms /u τ collapses for z + ≤ 15 and reaches its largest value for z + 15, while the v rms /u τ and w rms /u τ profiles do not collapse. However, they all agree with the analytic form imposed by the no-slip condition near the bottom boundary [START_REF] Pope | Turbulent flows[END_REF] u rms b 1 z, v rms b 2 z, w rms c 3 z 2 , (

where b 1 , b 2 , c 3 are constant coefficients. Fluctuations of the three velocity components are also significant within the logarithmic layer [START_REF] Lam | On the condition of streak formation in a bounded turbulent flow[END_REF][START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF], where they increase monotonically as the imposed shear increases. They remain proportional to the turbulent kinetic energy with a constant of proportionality that does not change with u τ .

A similar analysis can be done near the top surface (0 < z + s < 80). In that region (see Fig. 18), u rms /u τ profiles of cases ISC 0 , ISC 1 and ISC 3 reach a local peak at the surface, where momentum is injected by the imposed shear. This is an expected behavior, since a Taylor expansion of the various stresses near the top boundary implies that

u rms a 1 + c 1 z 2 , v rms a 2 + c 2 z 2 , w rms b 3 z, ( 2 2 )
where a 1 , a 2 , b 3 , b 3 and c 2 are constant coefficients. Similar to what was observed for the bottom boundary, fluctuations are also significant within the logarithmic layer. The rms fluctuations remain proportional to the turbulent kinetic energy as for the bottom wall, yet the constant of proportionality varies with u τ for the free surface.

D Structures

In this section, we present some qualitative pictures of the flow structure observed at horizontal planes located near the bottom wall (Fig. 19) and near the top surface (Fig. 20).

Fig. 1

 1 Fig.1Sketch of the flow configuration. The fluid domain is rectangular of length L x , width L y and height Lz . The origin is located at ( x, ŷ, ẑ) (0, 0, 0) so that bottom and surface are situated, respectively, at ẑ -ĥ

Fig. 2

 2 Fig.2Behavior of viscous and Reynolds stress (both normalized by the total mean shear stress) as a function of the wall-normal coordinate. Panel a: Situation near the bottom wall, as a function of z + . Panel b: situation near the surface, as a function of z + s . Case FSC is not computed at the surface, and therefore, it is plotted only in panel (a)

y 1 δFig. 3 Fig. 4

 134 Fig. 3 Fluctuation u of the instantaneous streamwise velocity, normalized by max x, y | u |, on the plane z + 5 located near the bottom boundary. Snapshots have been taken after reaching the statistically steady state. Only a portion of the numerical domain is shown to highlight the structures. The color bar is presented in quadratic scale. Top left figure, ISC 0 ; Top right figure, ISC 3 ; Bottom figure, FSC

  9 left) and (ii) the average Fourier spectrum of the vertical velocity component normalized by the surface shear velocity u τ s as a function of the wavelength in surface units λ + y, s 1 δ νs λ y at z + s 5 (see Fig. 9 right). One observes a dominant wavelength only for ISC 0 , ISC 1 and ISC 3 at λ y 0.45, 0.57, 0.48, respectively. This

Fig. 5 3 Fig. 6 30 Fig. 7

 536307 Fig. 5 Streamwise velocity u normalized by max y, z |u| in the transversal plane yz at the lower half of channel z + ∈ [0, 80]. Top panel, ISC 0 ; Central panel, FSC; Bottom panel, ISC 3

Fig. 8 Fig. 9

 89 Fig. 8 Instantaneous contour maps of the fluctuation u of the streamwise velocity, normalized by max x, y | u |, at the horizontal plane z + s 5, in surface units (but for FSC, in which the plot is located at 1-z δ ν 5). The color bar is presented in quadratic scale. Top panel, ISC 0 ; Central panel, FSC; Bottom panel, ISC 3 . Instantaneous snapshots are taken after a statistically steady state is reached

Fig. 10

 10 Fig. 10 Instantaneous contour maps of the horizontal velocity divergence ∇ H • u at the surface normalized by max x, y |∇ H • u| for the different cases. The color bar is displayed in quadratic scale. Top panel: ISC 0 ; Central panel, FSC; Bottom panel, ISC 3 . Instantaneous snapshots are taken after a statistically steady state is reached

  123

Fig. 11 Fig. 12

 1112 Fig. 11 Compressibility factor C (13) on the top half of the domain for all the simulated cases

Fig. 14 a 2 Fig. 16

 14216 Fig. 14 a Mean streamwise velocity profile u ; b mean shear stress τ ; and its decomposition in viscous (c) and Reynolds stresses (d). All cases are shown

Fig. 17 Fig. 18

 1718 Fig. 17 Root mean square of a streamwise, b spanwise, c wall-normal velocities, normalized by the friction velocity u τ and computed near the bottom boundary

Fig. 19 Fig. 20

 1920 Fig. 19 Fluctuation u of the instantaneous streamwise velocity, normalized by max x, y | u |, at the horizontal plane z + 5 located near the bottom boundary. Instantaneous fields are taken at one representative time inside the statistically steady regime for a ISC 0 ; b ISC 0 1; c FSC; d ISC 2 . e ISC 3 . Color bar is presented in quadratic scale

  Re. For a standard Poiseuille flow, the shear at the surface is opposite to that at the bottom τ s -τ w , and thus, u τ 1. Finally, a Reynolds number Re τ s u τ s Re, based on the friction velocity u τ s ≡ √ | τ s | at the surface.

	ρ | ˆ | ν Several other Reynolds numbers can be introduced. First, a bulk Reynolds number, .
	Re m 2u m Re,	with a bulk velocity u m ≡	1 2	1 -1	u dz,	( 8 )
	where angular brackets denote averaging in time and over the homogeneous directions x and y. Second, a friction Reynolds
	number Re τ ≡ ĥ ûτ /ν, based on the bottom friction velocity ûτ ≡ τw ρ . For FSC, the dimensionless bottom friction velocity is equal to u τ ≡ ûτ /[U ] √ 2 so that Re τ √ 2Re. For the ISC, u τ √ 2 + τ s with τ s ≡ τs ρ[U ] 2 leading to Re τ √ 2 + τ s

Table 1

 1 Mean flow parameters: the Reynolds number; the mean bulk velocity u m ; the associated bulk Reynolds number Re m ; the measured square of bottom friction velocity u 2 τ ; the associated friction Reynolds number Re τ ; the square of surface friction velocity u 2 τ s ; and the associated surface friction Reynolds number Re τ s

		Re	u m	Re m	u 2 τ	Re τ	u 2 τ s	Re τ s
	FSC	180	25.35	9.127 × 10 3	1.955	254.6	0	0
	ISC 0	180	13.80	4.968 × 10 3	0.843	165.3	1.157	193.6
	ISC 1	180	19.09	6.875 × 10 3	1.286	204.1	0.714	152.1
	ISC 2	180	27.24	9.807 × 10 3	2.249	269.9	0.272	93.8
	ISC 3	180	33.30	1.199 × 10 4	3.209	322.4	1.201	197.3

Table 1

 1

		ĥ(cm)	ν (cm 2 /s)	Ûm (cm/s)	ûτ (cm/s)	Ŵ / ĥ	Re m	Re τ
	Present (FSC)	3.35	9.5 × 10 -3	26.0	1.70	-	9127	594
	Zhong et al. [35]	3 . 3 5	9 .5 × 10 -3	29.2	1.81	9.0	10297	638

Table 3

 3 Summary of the main simulation parameters: the mean streamwise velocity at the surface V s ; the Reynolds number; the grid resolution in the streamwise direction in wall units x +x × Re τ , in the spanwise direction y + y × Re τ and in the normal wall direction z + z × Re τ (minimum and maximum values are indicated); the number of points along each direction, N x × N y × N z ; the time step t

		V s	Re	x +	y +	z +	Grid resolution	t +
	FSC	-	180	12.5	6.11	0.019 -3	256 × 256 × 257	1.5 10 -4
	ISC 0	0	180	8	3.96	0.01 -2	256 × 256 × 257	10 -4
	ISC 1	10	180	10	4.9	0.015 -2.5	256 × 256 × 257	10 -4
	ISC 2	35	180	13.2	6.5	0.02 -3.2	256 × 256 × 257	10 -4
	ISC 3	50	180	15.8	7.7	0.0024 -3.9	256 × 256 × 257	10 -4

Acknowledgements One of the authors (Jair Reyes) wish to acknowledge funding by a Mexican fellowship from Conacyt. We also wish to acknowledge Dr. Michael Berhanu for fruitful discussions.

Funding Open access funding provided by TU Wien (TUW).

Data Availability Statement

The authors provide in the manuscript and in the appendix all data necessary to understand, evaluate, replicate and build upon the reported research. Possible complementary data can be found in the cited references by the authors.


Antoneta Iuliana Bratcu

Iulian Munteanu
email: iulian.munteanu@gipsa-lab.grenoble-inp.fr

Olivier Sename
email: ‡olivier.sename@gipsa-lab.grenoble-inp.fr

Delphine Riu
email: delphine.riu@g2elab.grenoble-inp.fr

Hybridization of Energy Storage Units for Energy Transition Applications: A Focus on Control and Coordination Approaches

Keywords: transition, smart grid, control of microgrids, hybridization of energy storage units, robust control of multi-source electrical-energy systems

Hybridization of energy storage units is a topic of interest in the nowadays context of transition towards more sustainable ways of energy production and consumption. Extensive research continues to be devoted to this topic, whereas various technical solutions have already been successfully implemented for a plethora of applications. This paper aims in a first step at overviewing the necessity of storage units in an increasingly renewable-based, decentralized energy production, then at explaining the role of hybridization of different storage technologies for more versatility and flexibility in what is generally known as a smart grid context. Then, the focus is on how to make the different energy storage sources to optimally and robustly cooperate towards a common goal. It is about obtaining a "fusion" of different heterogeneous sources. The answer to this problem obviously requires advanced control approaches being employed. An overview of most effective and widely used control strategies is envisaged -accompanied by some illustrative application results -out of which robust control techniques are given a special attention. This paper ends by attempting a look towards the future, namely by identifying some open questions and research directions worthy to further investigate.

I. INTRODUCTION: CONTEXT OF THE GREEN ENERGY TRANSITION

A. Necessity of energy storage

The green energy transition is nowadays recognized as a necessity in response to growing environmental concerns resulted from depletion of fossil fuels. Practically all types of human activity have lastly benefited of advantages of renewable-energy conversion systems (RECS), enabled by the rapidly developing power conversion technology towards mitigating RECS unpredictability due to their stochastic dependence on natural conditions [START_REF] Ellabban | Renewable energy resources: Current status, future prospects and their enabling technology[END_REF], [START_REF] Blaabjerg | Overview of control and grid synchronization for distributed power generation systems[END_REF]. Indeed, RECS grid integration is challeging, mainly due to voltage and frequency instability, poor power quality and load following [START_REF] Kyritsis | Evolution of PV systems in Greece and review of applicable solutions for higher penetration levels[END_REF]. In this context, energy storage systems (ESS) come as a solution to increase RECS dispatchability, operation effectiveness and grid penetration [START_REF] Chong | Hybrid energy storage systems and control strategies for stand-alone renewable energy power systems[END_REF]- [START_REF] Hlal | Optimum battery depth of discharge for offgrid solar PV/battery system[END_REF]. Various ESS technologies include electrical, electro-chemical, chemical and mechanical storage systems, out of which batteries with different chemistries are now increasingly interesting, especially for grid stabilization purposes. An analysis of ESS and battery features, limitations and suitable applications is given in [START_REF] Babu | A comprehensive review of hybrid energy storage systems: converter topologies, control strategies and future prospects[END_REF].

The global energy storage market has continuously expanded in the last decade, partially as a result of national policies and subsidies provided by governments of developing nations to reach their eco-friendliness and sustainability goals. This is suggested by the evolution of their installed power and capacity (Fig. 1 [START_REF]Global Energy Storage Database[END_REF]), where the historically prevalent pumped hydro storage begins to leave more place to lithium-ion batteries. Indeed, the last ten-year cummulative sum of both installed rated power and capacity for the electro-chemical battery and chemical storage exhibits a steady exponential increase [START_REF]Global Energy Storage Database[END_REF]. The adoption of electro-chemical and chemical storages for various grid applications, as shown in Fig. 2 [START_REF]Global Energy Storage Database[END_REF], suggests again the growing interest of batteries, in particular redox flow and lithium-ion technologies.

B. Benefits of energy storage hybridization

Decision of integrating energy sources of different types and technologies must obviously be based on the conclusion that theses sources behave better together than each alone. Indeed, as ESS are limited either by their power or energy capacity, a combination of several technologies may be suitable with regard to employ each storage respectfully of its specificity [START_REF] Chen | Progress in electrical energy storage system: A critical review[END_REF]. A synthetic vision of HES adoption in energy grids is given in Fig. 3, emphasizing the HES great potential of improving grid flexibility and resilience [START_REF] Babu | A comprehensive review of hybrid energy storage systems: converter topologies, control strategies and future prospects[END_REF]. Thus, a key feature supporting the storage hybridization decision is the complementarity between the sources, whose meaning will be detailed next and which can be exploited as a first good premise, such as to lead further to their coordinated -possibly optimal with regard to some well-defined criterion -behaviour, whose guarantee can only be ensured by control action [START_REF] Dizqah | Standalone DC microgrids as complementarity dynamical systems: modeling and applications[END_REF]- [START_REF] Rajabzadeh | Dynamic modeling and nonlinear control of fuel cell vehicles with different hybrid power sources[END_REF].

The interconnection topology of HES within a grid/microgrid determines the controllability, dynamic Fig. 1. Installed rated power and capacity per year for all types of storage: electro-chemical battery and chemical, electro-mechanical and thermal [START_REF]Global Energy Storage Database[END_REF]. Fig. 2. Installed rated power of the electro-chemical battery and chemical storage per grid application [START_REF]Global Energy Storage Database[END_REF].

performance, efficiency and lifetime of the involved storages. In particular, connection by means of power-electronic converters ensures the degrees of freedom required for a flexible control and operation, unlike direct grid connection, which is in turn cost effective and relatively not complex. Fig. 4 [START_REF] Babu | A comprehensive review of hybrid energy storage systems: converter topologies, control strategies and future prospects[END_REF] presents an active topology -materializing a Fig. 3. A vision of the roles and benefits of HES within a grid [START_REF] Babu | A comprehensive review of hybrid energy storage systems: converter topologies, control strategies and future prospects[END_REF]. so-called common DC link that ensures the power transfer towards the AC grid/microgrid -where both HES and HPS (high-power storage) are fully controllable by means of dedicated DC-DC converters. This structure can be particularized as a semi-active topology if the HES converter is absent, or as a passive one if both DC-DC converters are absent (neither storage is controllable). Fig. 4. Active power-electronic interconnection topology in a multi-source system, giving by particularization either passive or semi-active topology [START_REF] Babu | A comprehensive review of hybrid energy storage systems: converter topologies, control strategies and future prospects[END_REF].

II. CONTROL AND COORDINATION WITHIN HYBRID-ENERGY-STORAGE (HES)-BASED SYSTEMS

A. Complementarity-based problem statement

To fix ideas, active-interconnection-based configurations of a generic multi-source microgrid are considered next, like the one in Fig. 5, which is equipped of classical -Diesel generator, thermal generator -as well as renewable sources -solar pannels, wind turbine -of storage units of various technologies -batteries, supercapacitors (ultracapacitors), fuel cells -and also of (possibly both) AC and DC loads. Partial versions of this configuration are also illustrative. Such a microgrid can operate in either grid-forming (stand-alone or autonomous) or grid-connected mode. In either case, sources' power flows are managed by a Power Management System (PMS), mainly consisting of the control algorithm ensuring proper coordination of sources towards a common goal. Power flow management acts at a shorter time scale than the Energy Management System (EMS); hence, PMS is placed nearer to the plant and lower in a possible multiple-layer control hierarchy. In particular, one of the PMS functions is to ensure the effective coordination between sources, such as to exploit their complementarity. Thus, the complementarity between two power source technologies will further specifically refer to their relative dynamic abilities in providing energy in a longer or a shorter term. Hence, complementarity supposes dynamic behaviours at different time scales; its dynamic substrate enables complementarity to justify and be also used for control and coordination purposes.

Dynamic specialization of storages is best emphasized by Ragone plot [START_REF] Kuperman | Battery/ultracapacitor hybrids for pulsed current loads: a review[END_REF] -see, for example, Fig. 6 for a quite complete view. Thus, dynamic specialization range of each storage can be identified through the relation between its energy density (specific energy) and its power density (specific power). A large span of storage technologies is illustrated: from high-specific-energy storage sources, such as fuel cells -able to provide energy in a long term and hence called "energy sources", so specialized in the low-frequency range of energy variations -to the high-specific-power sources, such as ultracapacitors (supercapacitors) -able to provide energy very quickly, hence called "power sources" and specialized in the high-frequency range. Different electrochemical storages (e.g., batteries) or mechanical storages (e.g., flywheels) are placed in between the two extremes. Notion of characteristic (own) frequency of a source/storage quantifies more precisely its dynamic specialization range:

f p = ρ p /ρ e , (1)
where ρ p is the power density, in W/kg and read on the abscissa of Ragone plot, and ρ e is the energy density, in J/kg and read on the ordinate. Characteristic frequency values for some of the main storage technologies are indicated in Fig. 7.

A given source behaves best when operated in its specialization range -i.e., when it must supply power variations placed around its own frequency -as its reliability and lifetime are best guaranteed this way. For sake of reliability, a high-energydensity source, also identified as main source can be suitably complemented -otherwise said, hybridized -by a high-powerdensity one, also named auxiliary source, in order to protect the first from possible harm due to unappropriate exploitation outside its specialization range (Fig. 8). Initial sizing and rating of a HES system rely upon two classes of indicators for hybridization, namely: power potential of hybridization (PPH), based on the possibility of reducing the main source's size -i.e., rated power, weight, cost, etc. -due to the presence of the auxiliary one, and energy potential of hybridization (EEH), based on taking into account variation spectrum of load power and also regularity of charging/discharging cycles of the main source if this is a storage one.

In a hybrid main-auxiliary source tandem the auxiliary source will supply the high-magnitude fast variations of power demand, which the main source is not able to supply without some risks for its longer-term operation. The simple decision to put together such two types of sources is not enough to guarantee their realiability-aware operation, a multi-source PMS being thus required to suitably coordinate them. Irrespective of the underlying method, this action relies upon dynamic or, equivalently, frequency separation of sources. Moreover, in order to address all optimization possibilities, PMS design may already be considered in the phase of initial system sizing, an approach that is called system-control co-design.

B. A brief look at the PMS design state of the art

Next, the focus will naturally be on hybrid (multi) energy storage used in microgrid applications. In a multi-source microgrid, design of EMS and PMS can be formalized as control problems, namely to ensure a desired load power sharing between the different sources, while preserving reliability of all sources and extending their remaining useful life by exploiting their complementarity [START_REF] Iannuzzi | Use of supercapacitors, fuel cells and electrochemical batteries for electric road vehicles: a control strategy[END_REF]. More precisely, the interest will be on the DC part of an active-interconnected microgrid (Fig. 4), i.e., the common DC link, which can easily accomodate fully controllable multiple storages, namely as current sources. This is a topic related to DC microgrids, which have lastly gained growing interest vs. AC ones, mainly because they do not need synchronization. Autonomous applications in avionic, automotive and marine industries, as well as power supply of remote areas exploit this advantage, among others [START_REF] Guerrero | Advanced control architectures for intelligent microgrids -Part I: Decentralized and hierarchical control[END_REF], [START_REF] Justo | AC-microgrids versus DC-microgrids with distributed energy resources: a review[END_REF].

A plethora of methods can be used to design efficient PMS. Model Predictive Control (MPC) -with its nonlinear and/or stochastic versions -may be employed in a multivariable control framework for stand-alone microgrids, with or without renewables, as it can take account of multiple constraints [START_REF] Dizqah | Standalone DC microgrids as complementarity dynamical systems: modeling and applications[END_REF], [START_REF] Schuler | L 1 -optimal control of large wind turbines[END_REF], [START_REF] Xie | An energy management strategy based on stochastic model predictive control for plug-in hybrid electric buses[END_REF]. However, need of an adequate plant mathematical model and of a priori knowledge of load variation to predict the behaviour appear as drawbacks. For example, in the case of electric vehicles (EVs), the resulted PMS depends on the driving cycle. More sophisticated approaches, such as nonlinear [START_REF] Alnejaili | Sliding mode control of a multi-source renewable power system[END_REF], [START_REF] Rajabzadeh | Dynamic modeling and nonlinear control of fuel cell vehicles with different hybrid power sources[END_REF] or multi-agent-based and games theory [START_REF] Khan | Multi-agent based distributed control architecture for microgrid energy management and optimization[END_REF] may also offer viable solutions.

Reliability-aware operation of each source may appropriately be embedded into a PMS design, for example, "translated" and/or adapted from its data sheet, leading to multicriteria techno-economical optimization. Thus, reliabilityaware PMS design usually implements optimization methods that integrate complex degradation and ageing models of storage units [START_REF] Martel | Optimal economybased battery degradation management dynamics for fuel-cell plug-in hybrid electric vehicles[END_REF]. The HES multi-time-scale nature justifies hierarchical control structures [START_REF] Xiandong | Hierarchical energy management system for multi-source multi-product microgrids[END_REF], composed of dynamically separated layers, according to sources being exploited in their "specialization" range of variations (see Ragone's taxonomy in Fig. 6), respectful of their reliability [START_REF] Castaings | Decomposed energy management of a multi-source fuel cell vehicle using energetic macroscopic representation[END_REF]. The idea is that highpower-density sources (e.g., ultracapacitors), to protect highenergy-density sources (e.g., fuel cells or batteries), whose current variations should be smoother. To this end, filtering ensures the most explicitly the dynamic separation, by restricting variations of current of each source within in a predefined frequency interval [START_REF] Silva | Application of a decoupling method based on online filtering technique for multisource electric vehicles[END_REF], [START_REF] Florescu | Adaptive frequency-separation-based energy management system for electric vehicles[END_REF], while an optimal LQG approach may also work, but limited to only two sources [START_REF] Florescu | LQG optimal control applied to on-board energy management system of allelectric vehicles[END_REF].

Continuing the idea of expressing the dynamic separation of sources by means of an optimization index based on their respective own frequencies leads to the generalization of the optimal approach for an arbitrary number of sources, namely by introducing more degrees of freedom. Thus, the PMS design by using weighting functions within an H ∞ robust control design allows the flexibility needed for suitably conditionning the sources' dynamic closed-loop behaviour, in addition to guaranteeing robust stability and performance. Possible combination with Linear Parameter Varying (LPV) techniques -e.g., in a polytopic approach that generalizes gain scheduling techniques -extends the domain of PMS effectiveness. The robust approach is detailed next.

III. ROBUST HES PMS DESIGN AND PERFORMANCE IN SOME ILLUSTRATIVE CASE STUDIES

A. Generic three-storage microgrid for PMS design statement Without loss of generality, a generic, three-storagefuel cell, battery and supercapacitor -microgrid operating as a power supply on board of an electric vehicle (EV) is here considered (Fig. 9 [START_REF] Nwesaty | Power sources coordination through multivariable LPV/H∞ control with application to multi-source electric vehicles[END_REF], [START_REF] Nwesaty | Robust energy management system for multi-source DC energy systems -real-time setup and validation[END_REF]). The storage active parallel configuration has the advantage of flexibility (sources' independent operation, facility to replace/add power sources, etc.). The number of sources, three, was chosen to illustrate the applicability of the approach to any number of sources (whereas only two sources could have been viewed as a particular case). The fuel cell acts here as the main source, it is unidirectional, so connected to a one-quadrant DC-DC converter. The two auxiliary sources -the battery and the supercapacitor, both bidirectional -must supply the relatively fast power demand variations and collect the reversed power (during braking phase), so they are connected each to a two-quadrant converter allowing charging/discharging. The electrical motor and its converter materializes the load. A key task of the PMS is to maximize benefits of hybridization, which suppposes that each source to operate in its dynamic specialization range.

First, a classical two-level control structure is adopted, where the upper-level LPV/H ∞ -based PMS provides references to the DC-DC converters' current control loops on a lower level (Fig. 10 [27]). The dynamic (frequency)-separation operation of the three storages is accompanied by some other constraints, so the following control objectives must be met:

• maintain the DC-link voltage constant, around a predefined V * dc setpoint, regardless of the load current variations. This goal applies to all configurations with different sources paralleled on a common DC link. It is this requirement that allows the HES PMS design statement as a disturbance-attenuation problem in the H ∞ framework;

• ensure dynamic separation of storages, i.e., each storage supplies power with respect to its own frequency according to Ragone's plot (see [START_REF] Ellabban | Renewable energy resources: Current status, future prospects and their enabling technology[END_REF] and Fig. 6). To this end, appropriate shaping of weighting functions associated to H ∞ control design is a crucial step in the design; • keep supercapacitor's state of charge (SoC) loosely around 50%, thus enabling continuously absorbing/providing power in response to instantaneous load power demand; • impose a desired steady-state power sharing between the "slow" storages, fuel cell and battery. In this way, the fuel cell may operate at a desired working point, e.g., the maximum-efficiency one. Steady-state behaviour could determine battery long-term charging cycle, according to its technology. Regarding the low-level control design, as the DC-link voltage V dc is kept constant at V * dc by the PMS, converters' currents have fairly linear dynamics, so classical PI tracking controllers can be used. The low level is designed to be fast enough to be neglected at the slower-dynamic upper level, so the following equalities hold next: I f c * =I f c , I bat * =I bat and I sc * =I sc . The three tracking low-level loops have the converters' duty cycles, α f c , α bat and α sc , as control inputs. In the LPV framework, these duty cycles compose the varying parameter vector:

ρ = [ρ 1 ρ 2 ρ 3 ] T ≡ [α f c α bat α sc ] T .
The sought-for LPV controller K(ρ) in Fig. 10 must ensure closed-loop stability, whereas minimizing the effect of load current, I load , acting as a disturbance on the DC-link voltage. Control design is thus tackled in the H ∞ framework applied to LPV systems, based on the so-called P-K form of control block diagram [START_REF] Lublin | H 2 (LQG) and H∞ control[END_REF] in Fig. 10, with the plant G(ρ) being extended with the weighting functions -denoted by Wwhose selection should obey the above-listed performance requirements. In Fig. 10 G

(ρ) plant's state vector is x = [V dc V 1 V 2 V 0 ] T (
with notations in Fig. 9), the exogenous vector, w, contains the disturbance signal, I load , and the reference signals, u is the control vector, y is the output vector composed of the errors in relation to references and z is the controlled output vector, containing the "shaped"by-weighting-function errors and control inputs. The detailed design procedure, including selection of weighting functions, can be found in [START_REF] Nwesaty | Power sources coordination through multivariable LPV/H∞ control with application to multi-source electric vehicles[END_REF].

Next, in order to provide a proof of concept of the robust HES PMS design methodology, real-time validation results are briefly presented in two illustrative case studies: a first one concerning an electromobility application and a secdond one in a frequency-regulation application in a stand-alone microgrid.

B. Case study 1: HILS validation of a PMS for a three-storage DC power supply on board of an EV Hardware-in-the-loop simulation (HILS) is nowadays widely known and used for control law rapid prototyping and preliminary validation, before any real-site implementation. Since its definition and introduction some thirty years ago for developing and testing control structures in the automotive industry using the MATLAB ® /Simulink ® software on a dSPACE TM real-time system [START_REF] Hanselmann | Hardware-in-the loop simulation as a standard qpproach for development, customization, and production test of ECU's[END_REF], HILS concept has increasingly been used and evolved to the status of a standard. This concept mainly consists in closing the loop around the controller embedded on a real device and the plant, which is only emulated and/or numerically simulated, because it is too expensive to use the real plant.

The first case study concerns a three-storage power supply on board of an EV and corresponds exactly to the configuration in Fig. 9. HILS validation results were obtained on a test bench designed in collaboration with FEMTO-ST Laboratory in Belfort, France, where real battery and supercapacitor, as well as a complete real-time fuel cell emulator developed at FEMTO-ST Laboratory [START_REF] Gao | PEM fuel cell stack modeling for real-time emulation in hardware-in-the-loop applications[END_REF] were used. For sake of reducing Fig. 9. Schematic of a three-storage -fuel-cell-battery-ultracapacitor -power supply system, detailing the equivalent electrical models of the storage units [START_REF] Nwesaty | Robust energy management system for multi-source DC energy systems -real-time setup and validation[END_REF].

Fig. 10. Closed-loop P-K form according to the LPV/H∞ formalism, where the extended plant contains the weighting functions [START_REF] Nwesaty | Robust energy management system for multi-source DC energy systems -real-time setup and validation[END_REF].

real-time implementation complexity, a reduced-order version of the controller was employed in real time, according to the order reduction method proposed in [START_REF] Nwesaty | Reduced-order LPV controller for coordination of power sources within multi-source energy systems[END_REF]. HILS validation required some adaptations, including in the selection of control references and weighting functions -see [START_REF] Nwesaty | Robust energy management system for multi-source DC energy systems -real-time setup and validation[END_REF] for details. Dynamic separation of storages is confirmed by time evolutions of the corresponding currents in Fig. 12.

C. Case study 2: HILS validation of a frequency-regulation PMS for a two-storage AC microgrid

Results presented within this subsection were obtained in a technology transfer project between 2018 and 2020 [START_REF]K-INF -fusion of electrical sources[END_REF], whose aim was to prove the practical feasibility of an H ∞based HES PMS solution and quantify its economical benefits. The use case was that of the primary frequency regulation in a two-storage-equipped stand-alone microgrid also containing both classical (Diesel generator) and renewable (PV) sources. The topology is a partial version of the one in Fig. 5, where a tandem battery-ultracapcitor materializes the HES.

Here, hybridization is justified by the large penetration rate of the renewables. The battery is employed to supply the load slow, long-term, low-frequency variations, while the supercapacitor, specialized in supplying fast, instantaneous, high-frequency (and possibly high-magnitude) load variations, is used as a "quick" backup storage that protects the battery. In this way, battery lifetime is saved, with some quantifiable economical benefits. Whereas the same HES PMS design methodology, generalizable to any number of storages/sources, was applied, however, for this case study some design and implementation details are purposely ommited, for reasons of intellectual property protection.

A test bench composed of a HILS-based demonstrator and a real battery pack together with its boost DC-DC converter -whose photo is presented in Fig. 13 -was built to allow submitting the battery pack to cycles relevant to the frequencyregulation use case, in order to further evaluate its aging parameters, out of which the residual capacity of discharge was the most relevant. To this end, the HES PMS version adapted to the frequency-regulation use case was embedded on a general-purpose TI C2000 microcontroller to allow predicting more than 50% increase of battery lifetime based on postcycling measurements, which suggests significant economical benefits of hybridization combined with robust HES PMSmore details can be found on the project website [START_REF]K-INF -fusion of electrical sources[END_REF].

The HES-PMS-managed microgrid closed-loop dynamic performance was challenged with the most unfavourable load power scenario, i.e., consisting of a series of step variations. Dynamic separation of battery and supercapacitor currents was concluded by analyzing the time evolutions of the concerned curents; explicit plots are skipped here for sake of space saving. On the other hand, a more general result was obtained in the frequency domain, which is shown in Fig. 14. It is about a histogram-based comparison of no-storage-vs. storage-based PMS: average of the microgrid frequency values are closer to the rated one (here 50 Hz) if the storage is present. This result offers a quantitative indication of the storage beneficial role in renewable-based microgrid control.

IV. CONCLUSION AND SOME PERSPECTIVES

The actual context of the green energy transition should rely upon extensive use of energy storage units, to support larger grid penetration and improve dispatchability of renewables. Hybrid energy storage (HES) has proven to improve grid ecofriendliness, flexibility and resilience. To this end, HES must be effectivey managed by Power Management Systems (PMS). This paper has argued the necessity of using advanced control methods for this purpose. Thus, details of systematically applying robust H ∞ control for the HES PMS design, as well as two illustrative case studies, have been presented.

To conclude, our proposed PMS design methodology requires input data representative of the dynamic specialization of storages, according to datasheet information and placement on the energy-density-vs.-power-density Ragone plot. A key notion is that of "own frequency" of a storage technology, as it precisely quantifyies its dynamic properties. Such an approach is coherent with a reliability-aware, thus sustainable design.

The HES PMS design has been stated as an S/KS mixedsensitivity problem in the H ∞ formalism, further allowing systematic robustness analysis. Analysis of the two HESbased microgrid case studies has shown that "fast" storages e.g., supercapacitors -are robustly well-performing in both improving overall closed-loop performance -in primary frequency regulation -and also in smoothing current variations of other "slower" storages -e.g., batteries -thus protecting them from over-sollicitation outside their "specialization" operation range, so prolonging their life and service time. In a more philosophical sense, the proposed PMS control design vision emphasizes the idea of harmonious fusion between different power sources/storages and can be regarded as an instance of a more general consensus principle.

Some interesting future research directions may be identified, as follows. Integrated into a microgird control, the proposed robust HES PMS design may require a computationally heavy LMI-based optimization problem being solved as the number of sources/storages increases. Even if scalability of the H ∞ -LPV design seems to be affordable for a small number of sources (two or three), it still needs some formal guarantees. Possible complexity alleviation by order reduction methods is a way to conforting the real-world application acceptability.

As the HES PMS must take account of roles of each storage within a given aplication, it must be compatible with some individual, possibly built-in, power magement systems -a Battery Management System (BMS) is a typical example, with the multiple functions it can provide in the case of some advanced versions, such as, for example, real-time estimation of state of charge and state of health. Thus, either clearer separation, or "fusion" of functions between the individual and global (HES) level may be envisaged. Also, from a system-control co-design viewpoint, specific economical benefits of storage service time extension by hybridization must be compared against initial investment and operation costs, to be able to really conclude on improving key performance indicators, such as CAPEX (capital expenditure) and OPEX (operation expenditure). Such an approach supposes complex, multi-objective optimization problems being appropriately formulated.

Fig. 5 .

 5 Fig. 5. General configuration of a microgrid with classical and renewable sources and various types of storage, coupled in parallel on a common DC link, supplying both AC and DC loads.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Ragone plot: dynamic specialization of different source/storage technologies.

Fig. 8 .

 8 Fig. 8. Simplified block diagram of a hybrid twostorage-unit supply system in a stand-alone microgrid.

  Fig. 11 exhibits an expanded view of the different parts composing the assembled test bench. The IFSTTAR (Institut Franc ¸ais des Sciences et Technologies des Transports, de l'Aménagement et des Réseaux -French Institute of Science and Technology for Transport, Development and Networks) driving cycle profile[START_REF] Nwesaty | Power sources coordination through multivariable LPV/H∞ control with application to multi-source electric vehicles[END_REF], representative of various urban driving conditions -comprising acceleration, deceleration, fixed speed and full brake -was chosen for being sufficiently rich in frequency content, such as to challenge the DC-bus voltage regulation and to illustrate the sources "fusion" in satisfying the power demand. The active load variation is controlled to correspond to the chosen driving cycle scenario by means of a rapid-prototyping platform based on the dSPACE TM MicroAutoBox ® II device.

Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .

 11121314 Fig. 11.Topology of the rapid-prototyping test bench at FEMTO-ST Laboratory in Belfort, France, used to validate the LPV/H∞ HES PMS[START_REF] Nwesaty | Robust energy management system for multi-source DC energy systems -real-time setup and validation[END_REF].

ACKNOWLEDGMENT

The authors would like to acknowledge work of their former Ph.D. students Adrian FLORESCU and Jean-Marc Waleed NWESATY, and of R&D engineers Sophie MIN-PICAULT and Julien MINET, who have contributed at HILS validations and proofs of concept. The authors also wish to thank their colleagues Alexandre RAVEY, David BOUQUAIN and Patrice LUBERDA, who are with FEMTO-ST Laboratory in Belfort, France, for their help with HILS validations.


Isabelle Viaud-Delmon
email: isabelle.viaud-delmon@ircam.fr

Olivier Warusfel

Peter Brugger

Marine Taffou

One participant of the Flight group could not complete the experiment. Auditory and body perception in microgravity

How is the conscious experience of our body impacted when information across sensory modalities are unusual, and how does it impact the localization of sound in space? Here we tested whether the absence of gravity induces an anomalous body experience and whether it interferes with spatial auditory perception, which strongly relies on a head-centered frame of reference. Participants were asked to indicate the lateralisation of auditory sources in different sensory conditions linked to vision and to gravity. The unique conditions of parabolic flight were used to ask whether gravity itself, or in combination with altered visual inputs, provides a critical cue for the emergence of distortions in body-centered frame of reference.

Methods

Auditory task

Participants were asked to indicate the lateralisation of auditory sources presented through headphones. Auditory stimuli were trains of gaussian white noise of 500 ms, pre-convolved with generic HRTFs at 20, 30, 40, -20, -30 and -40° in the azimuth and 0° in elevation in the frontal hemifield. Participants had to press a right or left key on a response pad to specify whether auditory stimuli were localized to the right or the left of their body.

Sensory conditions

The task was performed without or with vision, and with or without gravity for the Flight group (Fig. 2 and3). Four blocks were administered, each corresponding to a different sensory condition.

When visual information was available, participants were provided with an image of their body as if they would see themselves from an "out-of body" location, i.e. without the apparent left/right reversal inherent to a mirror image (Fig. 1A &B).

The four blocks were always presented in the same order: eyes closed in normogravity (EC-1G), eyes closed in microgravity (EC-0G), eyes opened visualising their inversed mirror image in normogravity (EO-1G), eyes open visualising their inversed image in microgravity (EO-0G). For the two 1G conditions, each block was completed during the steady phases of parabolic flight. For the two 0G conditions, the task was realised during the 20s weightless phases of parabolic flights.

Questionnaire measures

The sense of embodiment was studied through pre-and post questionnaires. At the end of the eyes open blocs (E0-1G and EO-0G), the Out of Body Experience (OBE) questionnaire (adapted from [START_REF] Botvinick | Rubber hands 'feel' touch that eyes see[END_REF] was administered to quantify the subjective experience associated with seeing one's self. Participants were asked to rate four different statements concerning their experience using a ten-point Likert scale.

Unusual perceptual experiences linked to the experiment were studied through a comparison of preand post scores of an adapted version of the Perceptual Aberration Scale (PAS) (Chapman, Chapman, and Rawlin, 1978), administered before and after the end of the entire protocol.

Participants

Twelve healthy human participants (7 female, mean age 35.00 ± 10.38) did the experiment during parabolic flights and constituted the Flight group. A Ground group of 35 healthy human participants (15 females, mean age 33.69 ± 9.19) did the whole protocol in normogravity. For this group, only two different sensory conditions were used since gravity did not vary. The two first blocks were done with the eyes closed in 1G (EC-1G), the two following blocs were done with eyes open visualising their inversed mirror image in 1G (EO-1G).

The scores of the questionnaires pre-and post-experiment assessing perceptual aberrations were significantly higher after the experiment compared to before in the Flight group. They were stable in the Ground group (Fig. 4).

Auditory task

The percentage of correct responses significantly differed according to the sensory condition in the Flight group, with the smallest percentage being found in the EO-0G condition (Fig. 6). They did not differ across conditions in the Ground group (Fig. 7).

Questionnaire measures

The OBE scores increased in the Flight group after the EO-0G bloc in comparison to the scores obtained after the EO-1G bloc, while they remained stable in the Ground group between the two EO-1G blocs (Fig. 5).

The sole manipulation of visual perspective linked to the body does not perturb the feeling of embodiment nor the frame of reference of auditory perception.

The frame of reference of auditory perception is influenced by gravity and visual perspective.

Figure 3 .

 3 Figure 3. Participant during the eyes open and normogravity condition (EO-1G), equipped with a head mounted display in which he is visualising himself.

Figure 2 .

 2 Figure 2. Participant equipped with a response pad and headphones in the eyes closed and microgravity condition (EC-0G).

Figure

  Figure 1A & B. In the eyes open conditions, a video camera captures the body of the participant from a distance of 150 cm, and this signal is sent to a software allowing to inverse the image before relaying it to the head mounted display worn by the participant.

FigureFigure 5 .

 5 Figure 4. Perceptual Aberration Scale scores in the Flight group (left) and in the Ground group before (PRE) and after (POST) the experiment.

Figure 6 .

 6 Figure 6. Percentage of correct responses in the Flight group in the four blocs (4 different sensory conditions).

Figure

  Figure 7. Percentage of correct responses in the Ground group in the four blocs (two sensory conditions).

Acknowledgements

This research was made possible by funding from the CNES. We are very grateful to the CADMOS and Sebastien Rouquette, to the CNES life science department and Guillemette Gauquelin-Koch, to the team of Novespace, to Pierre Denise and to our courageous participants.