Nonlinear hydrodynamics in a Mediterranean lagoon

. The paper addresses the application of the non-linear hydrodynamics model (RANS (Reynolds-averaged Navier–Stokes) equations) in a wide semi-enclosed Mediterranean lagoon (Berre lagoon), considering three natural forcing functions, i.e., sea tide propagating through a long narrow channel, wind and runoff. Main attention is focused to characteristic velocities (at free surface and bottom) and to free surface elevation associated to each of these three mechanisms, with special attention to the nearshore areas (i.e., in shallow water). The most interesting result concerns wind effects which, due to Berre lagoon bathymetry, give rise to downwind coastal jets, alongshore, in shallow water areas. Such coastal jets were never mentioned before in Berre la-goon literature.


Introduction
This paper concerns the hydrodynamics in a lagoonal ecosystem (called Berre lagoon) which was occupied, at the turn of the 20th century, by extensive Zostera noltii meadows (over 6000 ha), which shrinked to about 1.5 ha in 2004 due to environmental impacts of hydroelectric power and other anthropogenic developments as it is reported by Warner (2012).Over the last decades, there have been global declines in seagrass abundance in many places in the world, and management decisions aimed at protecting and restoring submerged aquatic vegetation have been taken in many places (Fonseca et al., 1996;Koch et al., 2006;Pickerell et al., 2005;Van der Heide et al., 2007;Cardoso et al., 2008;Koch et al., 2009, Shafer andBergstrom, 2010;Vacchi et al., 2012).Of course, it is known that there are a lot of abiotic and biotic factors, which can affect the losses of seagrasses and their lack of recovery.In particular, as kindly pointed out by one reviewer, changes in sediment properties, following the loss of seagrass, is probably an important factor to explain the lack of recovery.
In fact, the problem of the interaction hydrodynamicsmeadows needs a multiscale approach.In our case, the largest scale corresponds to the lagoon size itself (155 km 2 ): an intermediate would be the beach scale where meadow existed in the past (a few hundred square meters), and finally the meadow scale (a few square meters, for Berre lagoon) would have to be considered.But, the present research is limited to the largest scale, in which we will assume constant coefficients for both surface drag and bottom drag.As reported by Warner (2012), flow dynamics involving changes in water character, circulation and elevation in Berre lagoon, have not been studied in any detail up to now.So, the present paper is aimed to overcome this lack of knowledge and to better understand the impact of three main forcings (semidiurnal tide, wind and freshwater inflow) on the 3-D hydrodynamics in Berre lagoon: the circulation patterns and the bottom current in the nearshore areas (i.e., in shallow water).

Study area
Berre lagoon is one of the largest Mediterranean brackish lagoons (155 km 2 ; 0.98 × 10 9 m 3 ) in the South of France (Fig. 1).Its three main tributaries are situated in the northern part of the lagoon: an EDF (Electricité de France) industrial channel (exit of EDF hydropower station), with a runoff that can reach 250 m 3 s −1 during the peak exploitation period in winter, and two rivers (Arc and Touloubre), with mean runoff equal to 15 m 3 s −1 and 10 m 3 s −1 , respectively.Berre lagoon is connected to the Mediterranean Sea through a long and narrow channel, called Caronte.This complex Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.hydrosystem is subjected to semidiurnal tide, with amplitude up to 30 cm at the entrance of the channel.
The four control areas of the "Grand Etang" of Berre lagoon, marked by rectangles in Fig. 1, correspond to the places where benthic vegetation existed 50 yr ago and then disappeared completely (Figuerolles and Martigues, on the Western shore) or partially (Arc and "Pointe de Berre", on the Eastern shore).

Model statement
The numerical study is realized with MARS3D (3D hydrodynamic Model for Applications at Regional Scale) (Lazure and Dumas, 2008;Fiandrino et al., 2003) following the basic equations proposed by Blumberg and Mellor (1986).The software MARS3D is based on the system of incompressible Navier-Stokes equations in the classical Boussinesq approximation (ρ ∼ ρ 0 in momentum equation except for the buoyancy term) with the hydrostatic assumption: H L (Lazure and Dumas, 2008;Blumberg and Mellor, 1986).
Assuming that the free surface elevation is represented by the function z = η(x, y, t), and the bottom relief by z = −H (x, y), the model statement of MARS3D in coordinates (x, y, z) is where the pressure, P , is obtained by integrating hydrostatic equation ρg = −∂P /∂z from the bottom to the free surface; "prime" is a partial derivative u x ; u y ; u z = (∂u/∂x; ∂u/∂y; ∂u/∂z) .P (x, y, z, t) = P atm + ρ 0 gη + g 0 z ρ(x, y, z , t)dz , (4) V -velocity vector with horizontal components (u, v), f -Coriolis force, ν -coefficient of vertical turbulent exchange, ρ = ρ(S, T , P ) -water density, ρ 0 -reference density, µhorizontal eddy viscosity, F x = 2µu x x + µ u y + v x y , The model of MARS3D contains also two equations of thermodynamics for solving salinity and temperature fields: with F T ,S = k H (T , S) x x + k H (T , S) y y , where k H , k zcoefficients of horizontal and vertical diffusion.
The boundary conditions on the free surface z = η(x, y) are where (τ 0x , τ 0y ) -horizontal components the vector of the wind stress on the free surface, C p -constant, Q T -heat flux at the air-sea interface.
On the lower boundary; i.e., for z = −H (x, y), we assume that where (τ bx , τ by )-horizontal components vector of the tension stress on the bottom.
(τ bx , τ by ) = ρ 0 Cd B V (u, v) , where For the turbulence closure, the Prandtl model (Wilcox, 2004) was chosen to represent the vertical turbulent exchange coefficient.

Numerical modeling and mesh strategy
In MARS3D software, the water column is divided into layers in the transformed vertical "sigma" grid adapted to the bottom and free surface shapes: σ = (z − η)/(H + η); varying from −1 at the bottom −H , to 0 at the free surface η.For the discretization of Eqs. ( 1) -( 6) different grids are tested.A uniform grid for the horizontal axes and a nonuniform grid for the vertical axis are implemented as follows , where i, j, k -indexes for directions x, y, z; h x , h y , h σ -spatial steps; N x , N y , N σ -number of nodes along the two horizontal coordinate directions L x , L y , and along "sigma" direction, respectively.Concerning the horizontal directions, most of the numerical simulations are performed for h As a first variant of the vertical grid, the water column is simply divided into N σ −1 internal layers of nondimensional thickness h σ , and two limiting layers: BL (Bottom layer) and FSL (Free Surface Layer), of nondimensional thickness 0.5 h σ .The relationship between the thickness of the two limiting layers ε H = 0.5 H h σ and the number of node N σ is simply given as a function of the depth H in the water basin.For example, for the maximum depth (H = 9 m) the   Figure 4 shows the grid effect on the computed horizontal velocity components V = (u, v) in the water column, in the case of a flow driven by the wind only (oriented NNW, with a speed of 80 km h −1 ).The component u is positive in the sense of x increasing (from east to west), while the component v is positive in the sense of y increasing (from south to north).The results are given at a selected control place (point 1, in Fig. 3), where H = 9 m.
As a second grid variant, we implemented a refinement in both limiting layers, as indicated in Fig. 5 , in such a way that the nondimensional thickness of the two limiting layers is reduced by a factor 4, from 0.5h σ to 0.125h σ .In this variant, the number of vertical nodes is increased by 4: N σ = N σ +4.
To test simulation accuracy, computations are performed for N σ = 24 (i.e., for N σ = 20 and h σ = 1/20), for which the thickness of the two limiting layers is reduced to 5.625 cm at the most profound place (i.e., for H = 9 m).We concluded that this second grid variant, with N σ = 24, provides sufficient accuracy, and we decided to perform the present study with such a grid.

Main dynamical forcing mechanisms in Berre lagoon
In this paragraph, we perform a more detailed analysis of the effect of the three main dynamical forcings in Berre lagoon, i.e.: wind, tide, and freshwater inflow (runoffs).The following model conditions are considered: -sinusoidal tide (salt water inflow through a long and narrow channel, called Caronte, connecting the Mediterranean Sea to Berre lagoon), with amplitude up to 30 cm at the entrance of the channel; -constant wind velocity, for the two main orientations (NNW and SSE), ranging from 18 km h −1 to 100 km h −1 , according to the wind rose (Fig. 2).
We first analyze the impact of each forcing, separately, on the overall pattern of currents, and in the nearshore zones (typically 500 m× 500 m); in particular, the eastern nearshore zone called "Pointe de Berre" (shown in Fig. 1).In all cases, the initial condition used for the numerical simulation is the quiescent state (zero velocity everywhere) at t =0.For each set of meteorological conditions, it takes about two or three days for the numerical simulation to converge to a stable state, starting from the quiescent state.
Although our numerical simulation provides 3-D currents, only a few hydrodynamical characteristics: surface current, mean bottom current and barotropic currents (average along the vertical) will be shown and discussed in the present paper.To establish the harmonic behavior of the tidal effect in the whole lagoon, several tens of hours of numerical simulation starting from a quiescent state were necessary.The time evolution of the free surface elevation is presented in Fig. 6 for three characteristic points: (a) entrance of the Caronte channel (sea level), (b) exit of the Caronte channel into Berre lagoon, and (c) inside the Berre lagoon at the point 8 of Fig. 3 (close to the nearshore area of Pointe de Berre, ZB).The tide amplitude reaches about 7 cm at the exit of Caronte channel (red line) and about 6 cm near ZB (green line); that means that the tide inside decreases up to 5 times compared with the sea tide amplitude of 30 cm.We also observe a phase shifting -a time delay of tide propagation (free surface elevation maximum): inside the Berre lagoon of about 3 h at the exit of the Caronte channel (red line) and of about 4 h near the nearshore zone ZB (green line).
The present numerical results are well coherent with tidegauge observations reported by Ramade (1997) with tide amplitude decreasing up to 4-6 times, and a phase shifting of about 4 h.
In most part of the "Pointe de Berre" area, the velocity module is small enough (0-5 cm s

Wind effect in Berre lagoon
Due to the shallowness, currents and hydrodynamics of Berre lagoon are closely conditioned by the bottom topography, and wind affects the entire water column, as for many other Mediterranean lagoons (Perez-Ruzafa, 2011).Wind stress, which is caused by moving atmospheric disturbance, is known to have a major influence in lagoon water circulation.Because the scale of the lagoon body is small in comparison with the scale of cyclonic disturbances, the geostropic wind will be assumed to be uniform in the present numerical study; speed values of 80 km h −1 for each of the two main directions NNW and SSE will be considered; barotropic currents are shown in Figs.7 and 8, respectively.
It is observed that the (vertically integrated) barotropic current is maximal alongshore in the wind direction; it presents two principal gyres, longed downwind, which divide the "Grand Etang" in two parts, with opposite current to the wind direction in the deeper central part and two codirectional currents, stronger, along the eastern and western shores.
For a SSE wind of 80 km h −1 , Fig. 8 shows strong downwind currents along the eastern and western shores.Analogically to the NNW wind (Fig. 13), there are two main gyres, but circulating in the opposite sense.
Barotropic current fields, in Figs.7-8, give a global (vertically integrated) information about current structure in the whole lagoon, but cannot reflect complexity of the velocity distribution in the vertical direction.As an example of such complexity, we present in Fig. 9a, b, at the 8 control points (P1-P8) shown in Fig. 3, the vertical profiles of the two horizontal components of the velocity; the component u being positive in the sense of x increasing (from east to west), while the component v is positive in the sense of y increasing (from south to north).A change of the current direction is observed at a depth of about 1-2 m at the points 1, 2, 4, 5 and 8 (i.e., more inside), and of about 3-4 m at the points 3 and 6 (closer to shore).Point 7 is special, with a larger velocity; it corresponds to the zone of water exchange between Berre lagoon and the Caronte channel.
Figure 9c shows the profiles of the velocity modulus for 4 nearshore control points PB, PA, PM and PF.In the bottom layer, we determined an averaged velocity modulus (mean value at the three vertical mesh points in the limiting bottom layer (i.e., at h = H /40 from the bottom).Strong enough bottom currents exist at the center of these four nearshore control zones: 13.7 cm s −1 at PB, 16.2 cm s −1 at PA, 11.2 cm s −1 at PM and 10 cm s −1 at PF. Due to these high enough bottom velocities, which correspond to strong enough shear stress, we can expect an impact on the sediment composition, and so, on the erosion resistance of the sediment (Ahmad et al., 2011).Indeed, the living Zostera noltii are generally growing on mud-sand sediments; but when these plants disappeared in the past, due to fresh water stressor, the sediment composition was probably changed by a decreasing of mud concentration.
More detailed results of nearshore velocity fields at the free surface and at the bottom are given, as an example, in a square domain ZB (Pointe de Berre) of dimension 500 m × 500 m around the control point PB.The velocity fields for a NNW wind of 80 km h −1 are given in Fig. 10a for the free surface and in Fig. 10b for the sea bed.We can observe that starting from the control point PB (at the centre of ZB area) in direction of the shore, the free surface velocity diminishes monotonically when approaching the shoreline; but on the contrary, the bottom velocity reaches values substantially higher (Fig. 10b) in some parts of the shallow water.Typically, the bottom speed can be of the order of 0.2 m s −1 at PB (situated at 250 m from the shore), but becomes twice as fast at about 100 m of the shore.This observation is in agreement with the remark of Wiseman and Rouse (1980) about the importance of interaction with the topography (shallow bathymetry) mentionned by Allen (1975), in addition to the barotropic nature of the flow pointed out by Csanady (1971).

Influence of large freshwater runoffs
In this subchapter, the influence of large freshwater runoffs in the northern part of "Grand Etang" is analyzed in the case without wind or tidal effects.Two different values of EDF runoffs: 125 m 3 s −1 and 200 m 3 s −1 are considered.Additional inflows of the two rivers "Arc" and "Touloubre" (see Fig. 1), are also taken into account (15 m 3 s −1 and 10 m 3 s −1 , respectively).To analyze the impact of the EDF runoffs, we first consider the structure and intensity of the currents at the free surface layer for two selected runoffs (Figs. 11 and 12,respectively).
The flow coming from the EDF hydropower station is deflected by the western shore and divided in two branches; then the southern branch impinges the southern shore and is mainly deflected to the eastern direction.The contribution of river runoffs (Arc and Touloubre) is also visible.
In Figs.11 and 12, the global structure of the currents in the "Grand Etang" is the same, with two branches deflected by the western shore, but with stronger currents.In these cases, the deflection of the southern branch by the southern shore is stronger and gives rise to a counterclockwise gyre, with a significant returning current flowing along the eastern shore.
As in the previous subchapter, it is interesting to consider the vertical velocity profiles at some relevant places (8 control points in "Grand Etang", and 4 nearshore control points).downwelling current.It is to be noted that at the 4 nearshore control points, the current in the entire water column has the same direction than the surface flow shown in Fig. 12.In addition, we can see that the velocities induced by the EDF runoff taken separately (i.e., without wind or tide effect) in these 4 areas would not influence substantially bottom currents, except for the PB area (Fig. 13c).

Combined effects of strong wind and inflow
Of course, the three mechanisms, analyzed in the preceding subchapters, occur simultaneously and are independent of each other: the tides alternate approximately every 6 h; the wind is very often changing its direction and speed, while EDF runoff is seasonal.Now, on the basis on the results obtained in the previous subsections, we consider a situation in which two prevailing mechanisms (except tide) occur simultaneously and concurrently: SSE wind with a speed of 80 km h −1 and EDF runoff of 200 m 3 s −1 .For such a coupling, we could expect a larger velocity inside the eastern nearshore zones, like ZB.The influence of these two coupled forcings on the characteristic nearshore velocities (surface, bottom, average) is summarized in Table 1, which permits, in particular, a comparison between the results for SSE wind alone and coupled with large EDF runoff, respectively.We were expecting a co-current effect at PB. But, contrary to our expectation, the velocities are almost the same.That means that the co-current effect along the eastern shore is exactly balanced by the counter-current effect along the western shore.We observe almost the same balancing effect at the three other control areas, PA, PF and PM, for the same reasons.This means that the current is not simply the sum of the two contributions, confirming the strong nonlinearity of the system.

Discussion and conclusions
We have numerically investigated the nonlinear hydrodynamics induced by three main dynamical forcings (tide, wind and runoff) in Berre lagoon, a semi-closed hydrosystem which is subject to strong enough and frequent winds.This study was carried out as a contribution of a project of recolonization of protected seagrass (Zostera noltii) in some nearshore places in Berre lagoon.
The problem of recolonization of damaged seagrass is a very complex multiscale problem, which is the object of an always increasing number of studies, for both abiotic and biotic aspects, in many places in the world.Our contribution was concerning only one of the abiotic aspects (i.e., hydrodynamics) at the scale of the lagoon itself.That means that is was needed to make several assumptions; in particular, both free surface drag and bottom drag were assumed to be constant.
The use of the parallel version of MARS3D package permitted us to determine accurately the global influence of each mechanical forcing on the current structure in the whole lagoon and in the coastal areas.The circulation in the "Grand Etang" of Berre lagoon induced by the EDF inflow is a counterclockwise gyre.However, the pattern may be completely changed as a result of windforcing.For both main wind directions (NNW and SSE), the barotropic current is splitted into two separate gyres with coastal jets alongshore and windward.For strong NNW winds, the barotropic gyre is counterclockwise on the western side of the lagoon, and clockwise on the eastern side.For strong SSE winds, it is just the opposite.The main result, which was never mentioned before in Berre lagoon literature, is the existence of these coastal jets, and the fact that the bottom velocity is varying nonmonotonically in these shallow water areas when approaching the shore.Bottom velocity of the order of 40 cm s −1 , and faster, at a distance as small as 100 m from the shore, can be induced by a wind of the order of 80 km h −1 .
Such coastal jets with high enough bottom velocity, which correspond to strong enough shear stress, can be expected to have an impact on the sediment composition, and thus, on the erosion resistance of the sediment (Ahmad et al., 2011).Indeed, as the living Zostera noltii are generally growing on mud-sand sediments, we can expect that the sediment composition was probably changed by a decreasing of mud concentration.Such changes in sediment properties, following the loss of seagrass, could be an important factor to explain the lack of recovery.
It is known that the dependency of mud-sand composition on the erosion resistance of bottom sediments is a very open problem.So, it would be pertinent, and very challenging, in view of future recovery programs to study very carefully this question of hydrodynamics-sediment interaction, at least in the case of strong bottom current.

Figure 2 .
Figure 2. Distribution of the average wind over 10 min, station of Port de Bouc -10 m

Figure 4 . 6 Fig. 4 .
Figure 4. Profiles of horizontal velocity components V : (a) u and (b) v for different numbers 3 of layers σ N at point 1 in Fig. 3, for H =9m in the case of N-NW wind of 80 km/h, 24 4 hours after the beginning of simulation.5 6

Fig. 5 .
Fig. 5. Refinement schema for the vertical direction in the two limiting layers BL and FSL; red points -additional nodes in the two limiting layers.

Fig. 6 .
Fig. 6.Time evolution of the free surface elevation at three characteristic points: entrance of Caronte channel (sea level) -first line in black ; inside the Berre lagoon at the exit of Caronte channelsecond line in red ; inside the Berre lagoon at the point 8 of Fig. 4third line in green).

Figure 7 .
Figure 7. Barotropic currents and velocity modulus for N-NW wind of 80 km/h.

Figure 13 .Figure 13 .Figure 13 .
Figure 13.Vertical velocity profile for EDF runoff of 200 m 3 /s; u, v compon 7 control points inside "Grand Etang" (a,b), and (c) velocity modulus at the 4 nea 8 points PA, PB, PM and PF. 9 Figure 13 permits us to understand that the EDF runoff affects mainly the current near the free surface, for most of the control points in the central part of the Berre lagoon (3 m under the free surface) and the nearshore control points (1 m under the free surface).The surface velocity does not exceed 25 cm s −1 .At the points P1, P3 and P5 we can see a Nonlin.Processes Geophys., 20, 189-198, 2013 www.nonlin-processes-geophys.net/20/189/2013/

Table 1 .
Characteristic velocities, |V|surface, |V|bottom and |V|mean, at the four nearshore control points for different forcings.