An Unstructured Numerical Model to Study Wind-Driven Circulation Patterns in a Managed Coastal Mediterranean Wetland : The Vaccarès Lagoon System

The spatiotemporal structure of wind-driven circulation patterns and associated water exchanges can drive important bio-hydrodynamic interactions in shallow lagoons. The Vaccarès lagoon system is a complex shallow hydrosystem located in the central part of the Rhône Delta (France). It is internationally recognized as part of a biosphere reserve within the framework of UNESCO’s Man and Biosphere Programme, and as a RAMSAR site. Due to its frequent occurrence in this area, and considering the shallowness of the Vaccarès lagoon system, wind is assumed to play a major role in the hydrodynamic and biological processes. In this study, a hydrodynamic model was developed to investigate the structure of wind-driven circulations in the Vaccarès lagoon system, to provide insights into their role in transport and water exchange processes. The implementation and assessment (calibration and validation) of the model is presented first. Simulations were OPEN ACCESS


Introduction
In shallow lagoons, understanding the spatiotemporal structure of wind-driven circulation patterns and associated water exchanges is a crucial step towards assessing the bio-hydrodynamic interactions. The Vaccarès lagoon system, located in the Ile de Camargue, is an interesting system to study these wind-driven circulation patterns due to the shallow water levels and unusual complex geometry.
The Ile de Camargue basin is the central part of the Rhône Delta in the South of France, situated between the two branches of the Rhône River (see Figure 1). The aquatic ecosystem in the center of the Ile de Camargue is internationally recognized as a biosphere reserve within the framework of UNESCO's Man and Biosphere Programme, and as a RAMSAR site. The Ile de Camargue is of international importance for nesting, staging, and wintering water birds [1]. Among the various breeding species are Ardeidae (herons, bitterns, etc.), with very large numbers of Anatidae (ducks, geese, swans, etc.) occurring in winter.
The higher parts of this 800 km 2 area are devoted to agricultural land (420 km 2 ) consisting mainly of rice fields, whereas wetlands, brackish lagoons, and marshes occupy low-lying land [2]. This system includes three interconnected lagoons, which form the Vaccarès lagoon system. These interconnected lagoons consist of the Vaccarès lagoon, the Impériaux lagoon (IL), and the Lion/Dame lagoon (LDL) complex ( Figure 1). The Vaccarès lagoon is the only one that receives runoff from rice farming. Impériaux is the only lagoon exchanging water with the Mediterranean Sea, through the connection at La Fourcade (see Figure 1).
Tourism activities prevail in the western part of the delta. In the lowlands surrounding the lagoons, large freshwater marshes (25 km 2 ) are managed for waterfowl hunting, nature conservation, or exploited for their reed beds. Agricultural land borders the park to the north and southeast, and is mostly devoted to intensive flooded rice cultivation. The rice fields are irrigated during the crop period (from mid-April or early May, to September) by pumping stations taking water from both arms of the river Rhône. In the unpolderized agricultural area (112 km 2 ) in the northeastern part of the Ile de Camargue and around the Vaccarès lagoon, the runoff from the rice paddies is discharged directly into the lagoon system through two drainage channels (Fumemorte (FUM) and Roquemaure (ROQ), see Figure 1), whereas in the northern and southeastern parts, the runoff from 310 km 2 of polderized rice paddies is pumped back to the Rhône River or to the Mediterranean Sea [2]. In cases of very high precipitation, a portion of the runoff from the polderized rice paddies can be discharged directly into the lagoon system through the Rousty channel (ROU), see Figure 1. The volume of irrigation-induced drainage water discharged into the lagoons is about 480,000 m 3 per day during the crop period. This partially counterbalances the salinity increases in the lagoons that are due to high evaporation rates and seawater intrusion, but carries many anthropogenic substances into the protected area [3][4][5]. Thus, there is an ongoing debate on the management of water in this area. The Camargue natural regional park is separated from the Mediterranean Sea by a sea dyke. Sea-lagoon exchanges through the connection at La Fourcade (see Figure 1) are managed by 13 manual sluice gates of relatively moderate size (1 m). In autumn and winter, the gates are opened during northerly wind periods to decrease the water level in the lagoons. In spring, between one and three gates are opened to allow fish recruitment in the lagoons [6,7]. During sea storms, the sluice gates are closed. Although the seashore of the Rhône Delta is submitted to a micro-tidal regime (difference between lowest and highest value of astronomic tide is 0.42 m, see [8]), the "Vaccarès lagoon system" is not impacted by the tidal regime, due to the management of the sluice gates, mainly focusing on the decrease of the water level in the lagoons. The Vaccarès lagoon system is a complex system, exhibiting strong temporal and spatial variations in salinity, temperature, turbidity, oxygen, underwater irradiance, nutrients [9], and pesticide concentrations [4,5]. These spatial and temporal fluctuations in environmental conditions largely control the biogeochemical behavior of this system, and are the result of the interaction between freshwater and seawater inflows, rain, evaporation, and wind-driven forces, which vary over a wide range of time-scales and under particular geometry [10].
Experimental data indicate that strong prevailing winds can induce a change in water level of about 0.8 m considering two points in the northern part of the Vaccarès lagoon and in the southern part of the Impériaux lagoon [10]. Given that the average water depth is 1 m in the Vaccarès and about 0.5 m in the Impériaux and Lion/Dame lagoons, wind is assumed to play a major role in the circulation of water, salt, nutrients, contaminants and suspended sediments in the Vaccarès lagoon system. As evidenced in [11], in the Vaccarès lagoon, repeated lifting into suspension of sediment and horizontal transport over short time scales due to wind forcing are the major processes controlling the distribution of particles and adsorbed compounds. These in turn greatly affect the amount of light available in the shallow water column. Several studies have indicated that underwater irradiance in shallow ecosystems may vary significantly over different time scales as a result of wind forcing [12].
This study focuses on wind and has two main objectives: (i) to characterize the influence of wind on the circulation patterns in this complex hydrosystem; and (ii) to estimate the volume that is exchanged between the different lagoons in response to wind forcing. The model developed in this study is the first step towards enhancing our understanding of the biophysical interactions in the Vaccarès lagoon system (sediments, nutrients, pesticides, sea grasses, pathogen dynamics, etc.). This is the first study to explore the spatial structure of wind-driven circulation patterns for the whole Vaccarès lagoon system.
Given the small water depth values, it was assumed in a first approach that vertical movements are negligible, and that there is an absence of upwelling processes and stratification phenomena. The TELEMAC-2D two-dimensional vertically integrated model was then used [13,14]. This simplification has been adopted in the past by several authors [15][16][17][18][19]. The dependence of water density on temperature and salinity was not taken into account in this study.

Monitoring Data
During the last decade, an ongoing experimental dataset has been collected on the Vaccarès lagoon system, especially with regard to weather, water level, and water flow data.

Hydro-Meteorological Data
Wind (speed and direction), air temperature, and precipitation data were measured continuously at station B (see Figure 1). Wind (speed and direction) was also measured continuously at station A. All these data were averaged on a 60 min basis.
The Ile de Camargue basin is characterized by a Mediterranean climate, with mean annual precipitation (measured at meteorological station B, see Figure 1) of 615 mm (1963-2010) and a calculated (Penman method ) open water evaporation of 1477 mm (average value during the period  1988-2010).
An analysis carried out on wind data from stations A and B for the period 2002-2012 showed that the main wind regimes are north and northwest winds (directions 340° to 360°, see Figure 2a), with the average wind speed ranging from 4.9 to 7.6 m·s −1 . For these directions, see Figure 2b. Winds coming from the south (directions 120°-130°, see Figure 2a) are less frequent but can be intense during sea storms (average wind speed ranging from 6.5 to 7.9 m·s −1 , see Figure 2b).

Water Levels
Water levels were monitored continuously at 8 sites (see Figure 1), on a 15-min basis, since 2002 (float-operated Thalimedes Shaft Encoder with integral data logger from OTT Hydrometry). Three sites are located in the Vaccarès lagoon (1, 2 and 3, see Figure 1), two in the Impériaux lagoon (4 and 5, see Figure 1), two in the Lion/Dame lagoon (6 and 7, see Figure 1), and one in the Mediterranean Sea (8, see Figure 1). Water level gauges 1 and 2 are located respectively at the mouths of the Fumemorte and Rousty channels, see Figure 1. Data from the 8 locations are not always available in the same periods due to several problems (vandalism, thefts, devices temporarily out of order). (a)

Flow Data
At La Fourcade (lagoon-sea connection, see Figure 1), the discharge Q (m 3 ·s −1 ) through the gates is estimated using the following hydraulic equation: where Li is the width of the sluice number i (m), N is the number of opened sluices, H the experimental upstream water height above the sill (m), H' the experimental downstream water height above the sill (m), and Ki the discharge coefficient of the sluice number i (−). Values of the different Ki were determined with several flow measurement campaigns for each sluice.
At the outlet of the Fumemorte canal, the discharge was monitored continuously on a 30-min basis from 1993 to August 2008 using an automatic ultrasonic flow meter, which later stopped functioning. These datasets were used to develop and validate a rainfall-runoff model for the whole unpolderized drainage basin [10]. This model was used to calculate discharges from the unpolderized drainage basin into the lagoon system (see Figure 1).

Average wind speed
In addition, three gauging campaigns were conducted to measure water flows in several channels in the system: two campaigns in 2003 and 2010 in the channel between the Vaccarès and the Lion/Dame lagoons, and another in 2003 in one of the channels between the Vaccarès and the Impériaux lagoons (the most westerly), see Figure 1. The two gauging campaigns in 2003 were performed with an acoustic digital current meter (from OTT Hydrometry), and the campaign in 2010 with a surface mobile ADCP device, fixed on a small boat (OTT Hydrometry Qliner 2 ADCP Boat). The use of a surface ADCP in 2010 allowed us to perform faster measurements, in comparison with the acoustic digital current meter, limiting the risk of an inversion of the water flow in the channel during the measurement, due to a change in the wind direction.

Bathymetry
The bathymetric data were obtained from successive echo-sounding campaigns for the different lagoons, completed by one terrestrial LiDAR campaign for the outer edges of the lagoons and the inner islands. The total number of measured points was approximately 14,500.

General Overview of the Model
The two-dimensional mathematical model computing transient flows in the Ile de Camargue was constructed with TELEMAC-2D [13]. The TELEMAC-2D software solves the second-order partial differential equations for depth-averaged fluid flow derived from the full three-dimensional turbulent Navier-Stokes equations. This gives a system consisting of an equation for mass continuity and two force-momentum equations. The equations are averaged over the vertical by integrating from the bottom to the surface.
The averaged form of the continuity equation is: The averaged form of the momentum equations are: where u and v are the depth-averaged velocity components (m·s −1 ) in the x and y Cartesian directions; h is the water depth (m); Sh is the point injection or removal of fluid (m·s −1 ); νt is the coefficient of turbulent momentum diffusion (m 2 ·s −1 ); g is the gravitational acceleration (m·s −2 ); t is the time (s); Z is the free surface elevation (m); Sx and Sy (m·s −2 ) are the source terms of the momentum equation in u and v, respectively, including the Coriolis force, bottom friction, surface wind shear, and source/sink of momentum in the domain. In our application, TELEMAC-2D solves Equations (2)-(4) in two steps: (i) solution of advection terms; and (ii) solution of the propagation, diffusion, and source terms. The method of characteristics has been applied to solve the advection of u and v. The conservative PSI-scheme (Positive Streamwise Invariant scheme) has been applied to solve the advection of h. The propagation, diffusion, and source terms have been solved by the finite element method. Time discretization is used to eliminate non-linear terms from the equations. Variational formulation and discretization in space and time then transform the continuous equations into a discrete linear system whose unknowns are the values of u, v, and h at the nodes of the mesh. This system is then resolved by the GMRES method (Generalized Minimum RESidual method) [20,21]. Smagorinsky's formulation [22] has been used to parameterize the horizontal eddy viscosity.
Due to a lack of experimental data, spatial and temporal variations of atmospheric pressure were not taken into account in the model. Coriolis force was considered in the model as a constant source term of the momentum equations.

Mesh and Altitude
Equations (2)-(4) are solved on an unstructured grid, made with triangular elements. Thanks to the nature of this mesh, elements of various sizes can be fitted, thereby allowing high resolution in areas with small channels, and low resolution in areas where detail is not required.
The unstructured mesh generated for the Vaccarès lagoon system consisted of some 35,200 triangular elements of varying size and 19,700 nodes (see Figure 3). Special attention was devoted to the description of the different channels between the lagoons, and to the locations of the water level probes (see Figure 1). The triangular elements ranged in size from 300 m in the different lagoons to 1.5 m in water level probe areas. In the different channels between IL and LDL lagoon subunits, the size of the triangular elements was about 2 m, as illustrated in the Supporting Information.
The different lagoons of the Vaccarès hydrosystem at a water level of 0 m NGF are characterized by an average water depth of 0.5 m in the southern lagoons, 1 m in the Vaccarès (representing approximately 80% of the total water volume in the lagoon complex), with a maximum depth of 2.2 m. The French topographic reference (zero mNGF) was fixed in Marseille at the end of the 19th century by averaging local sea level records between 1885 and 1897 [23,24]. For a water level in the Vaccarès lagoon system of 0 mNGF, the volume of water is about 101 Mm 3 , and the area covered is about 110 km 2 .

Bed Friction
Bed friction is estimated using the Strickler's law, see Equations (5) and (6).
where α is the steepest slope of the bottom at the considered point, and K is the Strickler coefficient (m 1/3 ·s −1 ). The state of the sea grass meadow in the Vaccarès lagoon has been monitored annually since 1996 [25,26]. Regarding the characteristics of the bottom (type of sediment, sediment grain size, etc.) sediment samples were collected and analyzed in 2003 on 23 stations equally distributed in the Vaccarès [25]. However, very little data exists for the Impériaux and Lion/Dame lagoons. Due to this lack of information, it was decided to use a constant Strickler coefficient to model the bottom friction. This simplified approach must be improved in the future, with respect to the high spatial variability in the sea grass bed in this hydrosystem [25]. This Strickler coefficient was used as a calibration parameter of the water depth and the circulation between lagoons (final value of 40 m 1/3 ·s −1 ).

Boundary and Initial Conditions
Three liquid boundaries were considered to represent (i) the connection between the Mediterranean Sea and the Impériaux lagoon (IL, see Figure 1) and (ii) the two drainage channels (Fumemorte (FUM) and Roquemaure (ROQ), see Figure 1).
For the Fourcade liquid boundary, the imposed water flow was estimated at each time step using Equation (1), which has been implemented in TELEMAC-2D. Existing experimental data were used for the sea water level variations (see level gauge 8 in Figure 1) and the number of open sluices. The water level in the Impériaux lagoon is the water level calculated by TELEMAC in a pre-defined node close to the boundary. At the Fumemorte and Roquemaure boundaries, the water flows estimated by the rainfall-runoff model [10] were imposed. Other boundaries are rigid or slip boundaries.
Initial water levels were estimated using experimental data from the level gauges one to seven, see Figure 1. These seven point measurements were extrapolated to the different lagoons to obtain the initial water level conditions.

Rainfall and Evaporation
Daily potential evaporation was estimated using the Penman method [27]. In a first approach, this daily evaporation was then divided by 24 to obtain an approximation of hourly evaporation. Hourly precipitation data is measured directly at meteorological station B. At each time step, imposed evaporation and rainfall were estimated by interpolation of these hourly data.
where Uwind and Vwind are the components in the x and y Cartesian directions of the wind velocity (m·s −1 ), Cd is the drag coefficient (−), ρair the air density (kg·m −3 ), and ρ the water density (kg·m −3 ).
To account for the differences between winds recorded at the two meteorological stations A and B, the spatial interpolation described in Equations (9)-(12) is used: where UA, VA, UB and VB are the velocity components (m·s −1 ) of the wind velocity measured at stations A and B, respectively. XA, YA, XB and YB are the horizontal space coordinates of the two meteorological stations A and B (Figure 1), DA (m) the distance from the meteorological station A (XA,YA) to a computational node of coordinates (X,Y), and DB (m) the distance from the meteorological station B (XB,YB) to a computational node of coordinates (X,Y). The drag coefficient does not have a constant value. The following empirical relation was used to determine its value, adapted from [28]: user-specific wind speeds (m·s −1 ), used as calibration coefficients.

Areas Alternately Dry and Wet
Areas alternately dry and wet during a simulation were considered by using a specific computation of the free surface gradient as soon as the free surface intersects with the bottom of an element; see [20] for more details.

Calibration
The model was calibrated using experimental data from a 21-day period in May-June 2010 (Run 1: 24 May 2010 12 a.m.-14 June 2010 12 p.m.). During this process, adjustments were made to the bottom friction and Cd (see Equation (13)) coefficients to match the measured water levels in the lagoons and discharges into channel C1 (see Figure 1).  (13)).
Water level measurements from locations 1-4, 6 and 7 (see Figure 1) are plotted in Figure 4 together with the simulation results (no data were available for location five for this period). Initial water levels vary from 0.08 m NGF for location four to 0.15 m NGF for location six. The accumulated precipitation and evaporation were 6.6 mm and 143.6 mm, respectively. The performance of the model was analyzed using the four indices: mean absolute error (MAE), root-mean-square error (RMSE), normalized root-mean-square error (NRMSE) normalized as a percentage value of the average local water depth, and index of agreement (d) [29,30]. The definitions of these four indices are given in the Supporting Information.
For Run 1, the mean absolute error MAE, the root-mean-square error RMSE, the normalized root-mean-square error NRMSE and the index of agreement (d) are given in Table 1. There is reasonably good agreement between the observed and simulated water levels at all the locations, with the MAE varying from 0.01 to 0.03 m, the RMSE varying from 0.02 to 0.04 m, the NRMSE varying from 2% to 18%, and the index of agreement varying from 0.9 to 0.99. The model performs reasonably well during calibration.
Values of MAE, RMSE and d for the three lagoons (Vaccarès lagoon, Impériaux lagoon, and the complex Lion/Dame lagoon) of the Vaccarès lagoon system do not show significant differences, except for location 3 with the highest MAE, RMSE and NRMSE.
One series of eight flow measurements in channel C1 (see Figure 1) was performed on 27 May (surface ADCP device, see Part 2.1.3.). Water discharge measurements from this channel are plotted in Figure 5 together with the simulation results, after calibration. No experimental discharge data were available in the other channels for this period.
There is good agreement between the observed and simulated water discharges for channel C1, with a MAE of 0.43 m 3 ·s −1 , a RMSE of 0.66 m 3 ·s −1 , and an index of agreement of 0.94. It would be interesting to perform additional flow measurements in this channel to confirm these results.

Validation
After calibration, the model was validated by comparing simulations with water level data and, when available, water discharge data from October and December 2003.  Figure 1) are plotted in Figure 6 together with the simulation results (no data were available for location seven for this period). Initial water levels vary from −0.16 m NGF for location 1 to 0.02 m NGF for location four. The accumulated precipitation and evaporation were 1.6 mm and 39 mm, respectively.  For Run 2, the mean absolute error MAE, the root-mean-square error RMSE, the normalized root-mean-square error NRMSE and the index of agreement are given in Table 2.
There is reasonably good agreement between the observed and simulated water levels at all the locations (see Figure 6), with the MAE varying from 0.01 to 0.03 m, the RMSE varying from 0.02 to 0.05 m, the NRMSE varying from 3% to 34% and the index of agreement varying from 0.86 to 0.95. Values of MAE, RMSE and RMSE are less good for the locations four and five for which the model clearly overestimates the influence of wind. Differences of about 0.2 m can thus be observed once at location five and four times at location four, see Figure 6. Regarding discharges, no experimental data were available in channels C1, C2, C3, and C4 for this period.  Figure 1) are plotted in Figure 7 together with the simulation results (no data were available for location seven for this period). Initial water levels vary from 0.29 m NGF for the location four to 0.47 m NGF for the location three. The accumulated precipitation and evaporation were 6.4 mm and 7.6 mm respectively.
For Run 3, the mean absolute error MAE, the root-mean-square error RMSE, the normalized root-mean-square error NRMSE and the index of agreement are given in Table 3.
The observed and simulated water levels at all the locations agree well again (see Figure 7), with the MAE varying from 0.01 to 0.02 m, the RMSE varying from 0.02 to 0.03 m, the NRMSE varying from 3% to 5% and the index of agreement varying from 0.87 to 0.96. Values of MAE, RMSE and d for the three lagoons of the Vaccarès lagoon system do not show significant differences.
One series of flow measurements in channels C1 and C3 (see Figure 1) was performed on 15 December (acoustic digital current meter device, see part 2.1.3.). No experimental discharge data were available in channels C2 and C4 for this period. The two values of water discharge measured in channels C1 and C3 are 16.4 m 3 ·s −1 and 19 m 3 ·s −1 , respectively, for simulated water discharges of 16.1 m 3 ·s −1 in channel C1 and 17.2 m 3 ·s −1 in channel C3. There is rather good agreement between the observed and simulated water discharges for channels C1 and C3. The measured discharge in C1 (16.4 m 3 ·s −1 ) is about eight times higher than that used for the calibration (about −2 to 2 m 3 ·s −1 , see Figure 5).

Steady State Circulation Patterns and Water Levels
For real unsteady wind conditions, the circulation patterns continually change with the wind speed and direction, and a steady state is rarely attained in complex systems like the Vaccarès lagoon system. It is; however, instructive in a first approach to study and discuss results from hypothetical scenarios where the wind forcing is constant [17]. Two simulations were therefore performed for constant northwesterly and southeasterly winds. All system boundaries were specified as closed for these two simulations (no connection with the Mediterranean Sea nor with the unpolderized agricultural area, and no rain and evaporation). The initial simulated water level was 0.06 mNGF, corresponding to the mean water level in the Vaccarès lagoon system from 1993 to 2009 (unpublished data). Referring to the analysis carried out on wind data from station B for the period 2002-2012, the two most typical forcing winds simulated were (i) wind direction 350° and average wind speed 7.59 m·s −1 , and (ii) wind direction 130° and average wind speed 6.45 m·s −1 .
The steady state circulation patterns and water levels in the Vaccarès lagoon system for the most typical and constant northwesterly and southeasterly winds are shown in Figure 8a (direction of 350° and speed of 7.59 m·s −1 ) and Figure 8b (direction of 130° and speed of 6.45 m·s −1 ).
Under steady wind forcing, and for an initial water level of 0.06 mNGF, the circulation patterns and water levels in the Vaccarès lagoon system take about six days to reach an approximate steady state for the two tested wind conditions. This time scale is clearly longer than that of wind speed and/or direction changes. This means that in the Vaccarès lagoon system the circulation patterns are in a transient state and a steady state is rarely attained, except for an absence of wind over a period of several days.
For the two tested wind conditions, the water levels increase monotonically in the downwind direction across the different lagoons. Due to its geometrical characteristics, the Lion/Dame lagoon is the most impacted in terms of water level changes for the two wind conditions tested (see Figure 8a,b).
For northwesterly or southeasterly winds, depth-averaged flow velocities are relatively high and in the same direction as the wind along the shallows of the Vaccarès lagoon. In these areas, the surface wind stress and bottom friction tend to dominate hydrostatic pressure gradients. In deeper parts of this lagoon, there are typically slower counter-flows that oppose the wind-direction. In these areas, the axial hydrostatic pressure gradient can dominate the wind stress. These different flows along the shallows of the Vaccarès lagoon and in deeper areas are interlinked by circulatory gyres. In this lagoon, three main circulatory gyres are observed for the northwesterly wind (Figure 8a). These results are similar to those obtained by [11] in the Vaccarès lagoon for a northwesterly wind (direction of 315° and speed of 8 m·s −1 ). For the southeasterly wind, one circulatory gyre (anti-cyclonic rotation) is observed in the southern part of the lagoon, whereas a very large area of recirculation (cyclonic rotation) is observed in the northern part. In the southern part of the Vaccarès lagoon, there is an area with a bathymetric gradient (see Figure 3). Within a distance of about 1000 m, bathymetry changes from about −1.80 mNGF to about −0.6 mNGF. For both wind conditions, this change induces a pronounced increase in the velocities of the "north-to-south" currents and a pronounced decrease for the "south-to-north" currents (Figure 8a,b).
In the Impériaux and Lion/Dame lagoons, circulations are more complex. Some of the observed circulatory gyres and currents induced by the two simulated winds are represented in a very simplified approach. This is not an exhaustive representation of all the complex phenomena, gyres, and currents in this area, but it allows certain information to be deduced from these simulations.  Figure 8. (a) Steady state circulations and water levels in the Vaccarès lagoon system due to a northwesterly wind (direction of 350° and speed of 7.59 m·s −1 ); (b) Steady state circulations and water levels in the Vaccarès lagoon system due to a southeasterly wind (direction of 130° and speed of 6.45 m·s −1 ). Direction is defined as the direction from which the wind is blowing. Circulatory gyres are shown in a schematic and simplified way, with crosses for anti-cyclonic rotation and circles for cyclonic rotation. Only the main gyres are represented. For the zooms 1-4, a finer spatial interpolation of the current field was performed, and the scale used for the vectors is twice as high compared to the rest of the figure. Areas in black are no water areas.
In the Impériaux lagoon, for the northwesterly wind, two north-south currents can be identified. The first takes place along the western shore (western part of the island I1) where, generally speaking, water flows towards the south, with counter flows towards the north along the east coast of this island and the west coast of the island I2 (see zoom 1). The second current is observed in the area between the islands I2 and I3: water flows overall towards the north in the central deeper area between I2 and I3, with counter flows towards the south in the higher parts of this area (situated along the shores of these two islands (see zoom 2). In the northern part of this lagoon, several circulatory gyres are observed, associated with complex currents. Five of these gyres are shown schematically in Figure 8a.
For the southeasterly wind, a current can be identified in the northern part of the Impériaux lagoon, where water flows overall towards the west along the north coast of this lagoon and along the shore of the island I3. Interestingly, the two north-south currents observed for the northwesterly wind (see Figure 8a), are again observed with inversed orientations. The first takes place along the western shore (western part of the island I1), where water flows overall towards the north, with low counter flows towards the south along the east coast of this island (see zoom 3). The second current is observed in the area between the islands I2 and I3: water flows overall towards the south in the central deeper area between I2 and I3, with low counter flows towards the north in the higher parts of this area (mainly situated along the western shore of I3, see zoom 4).
In the Lion/Dame lagoon, for the northwesterly wind, the water flows overall towards the south along the western shore (western part of the island I7), and returns to the north along the east coast of this island. These flows are interlinked in the northern part of the lagoon by two circulatory gyres (anti-cyclonic and cyclonic rotations), see Figure 8a. For the southeasterly wind, the water flows overall towards the south along the eastern shore (eastern part of the island I7), and returns to the north along the west coast of this island. These flows are interlinked in the northern part of the lagoon by three circulatory gyres, as shown in Figure 8b. For the two wind conditions (Figure 8a,b) a circulatory gyre with anti-cyclonic rotation is observed in the southeastern part of the Lion/Dame lagoon.
Several smaller currents and gyres are induced by the two tested winds in the Vaccarès lagoon system, but are not represented in Figure 8a,b for the sake of legibility of this figure.

Unsteady Circulation Patterns/Water Levels Variations and Volume Exchanges between the Different Lagoons
To discuss the water levels variations in terms of the wind variability and to obtain information on the associated volume exchanges between the Vaccarès lagoon, the Impériaux lagoon, and the complex Lion/Dame lagoon, the real wind conditions of the period 24 May 2010 12 a.m.-14 June 2010 12 p.m. (the same period as for Run 1) were simulated. The initial simulated water level was 0.13 mNGF, corresponding to the mean water level measured in the Vaccarès lagoon system. All system boundaries were specified as closed for this simulation. For an initial water level of 0.13 mNGF, water volumes in the Vaccarès, the Impériaux and the "Lion/Dame" lagoons are respectively 122.1 Mm 3 , 31.6 Mm 3 and 6.9 Mm 3 .
Regarding the wind during this period, seven characteristics periods can be defined, as indicated in Figure 9. The first period ("P1" in Figure 9), from 24 May to 25 May, is characterized by southeasterly and southern winds, with a moderate wind speed (wind speed less than 5 m·s −1 , average value of 2.3 m·s −1 ).
The second period ("P2" in Figure 9), from 25 May to 27 May, is characterized by an alternation between southeasterly and northeasterly winds (direction between 0° and 180°). For this period, the average wind speed is 5.6 m·s −1 .
The third period (P3 in Figure 9), from 27 May to 29 May, presents a high amplitude in the wind direction (wind direction ranging from 0° to 360°). For this period, the wind speed is mostly less than 5 m·s −1 , with an average value of 3.7 m·s −1 .
The following period, from 29 May to the end of 3 June (P4a and P4b in Figure 9, identified by the first gray box), is characterized by northwesterly winds. For this period, the wind speed is mostly faster than 5 m·s −1 . Two sub-periods can be identified in period 4, the first one (P4a in Figure 9), from 29 May to 31 May with a wind direction close to 280°, and the second one, from 31 May to the end of 3 June (P4b in Figure 9), with wind directions close to 330°. Period 5, from the end of 3 June to the end of 7 June (P5 in Figure 9), is characterized by many changes in the wind direction (ranging from 0° to 360°). This is a transitional period linking period P4b (north-westerly winds) to period P6 (from the end of 7 June to 11 June), characterized by southeasterly winds (direction from 90° to 180°). For Period 6 (identified by the second gray box in Figure 9), the wind speed is mostly faster than 5 m·s −1 , with an average value of 8.3 m·s −1 .
The last period (P7 in Figure 9), which begins on 11 June, is characterized by many changes in the wind direction (ranging from 0° to 360°). The wind speed is mostly less than 5 m·s −1 , with an average value of 2.2 m·s −1 .
Several pieces of information can be derived from the Figures 10 and 11, which indicate for periods P1 to P7 the volume exchanges between the Vaccarès lagoon, the Impériaux lagoon, and the complex Lion/Dame lagoon (Figure 10), and the mean water levels variations in these lagoons ( Figure 11).
First, Figure 10 shows that the variations in the Vaccarès and the Impériaux lagoons are closely linked, with the volume in one increasing as the volume in the other decreases (Figure 10b,c). This relationship cannot be observed in the volume variations between the Vaccarès and the Lion/Dame lagoons (Figure 10d).
Due to the overall geometry of the Vaccarès lagoon system, water exchanges induced by wind between the different sub-lagoons mainly occur between the Vaccarès and the Impériaux lagoons. This is illustrated in Figure 10e, showing the exchanges between the Vaccarès lagoon and the Impériaux and Lion/Dame lagoons respectively. Regarding the Lion/Dame lagoon, wind-induced water exchanges mainly occur with the Vaccarès lagoon (results not shown in this paper). This is explained by the small size of the channel C4 compared to channel C1 (see Figure 3). These important results are expected to have an impact on the biogeochemical dynamics of the entire Vaccarès lagoon system. Another point is that the water levels variations in the Vaccarès, the Impériaux, and the complex Lion/Dame lagoons have different responses to wind variations, as illustrated in Figure 11.
Globally, for periods P1 and P2, there is a decrease in the water levels in the Lion/Dame lagoon. This increase is more pronounced for period P2, due to the higher wind speed (average wind speed is 5.6 m·s −1 for P2, and 2.3 m·s −1 for P1). In addition there is globally a decrease in the water levels in the southeast of the Vaccarès lagoon and of the Impériaux lagoon. Due to the small wind speed, the mean variations observed are not high (range from −0.04 to 0.02 m).
For period P3, there is globally a decrease in the water levels in the western part of the lagoons, and an increase in the eastern parts. Although the average wind speed is smaller than for the period P2 (5.6 m·s −1 for P2, and 3.4 m·s −1 for P3), the amplitude of the mean variation in the eastern part of the Lion/Dame lagoon is higher (0.04 to 0.06 m, in comparison with −0.04 to −0.02 m).
Globally, for the period P4a, with wind directions close to 280°, as for period P3, there is a decrease in the water level in the western part of the lagoons, and an increase in the eastern parts. However, due to the higher wind speed (average value 7.7 m·s −1 ), the mean variations are higher (see Figure 11, P4a). The highest decreases are observed in the Imperiaux lagoon, whereas the highest increases are observed in the Lion/Dame lagoon. For the 21 days simulated, the maximum volume in the Vaccarès lagoon, and the minimum volume in the Impériaux lagoon are observed for the period P4a. For this period, water flows in channels C2 and C3 are both from the Impériaux lagoon into the Vaccarès lagoon.

P3 P4a
P4b P5 P6 P7 Figure 11. Mean modeled water levels variations (m) in the Vaccarès, the Impériaux, and the complex Lion/Dame lagoons for periods P1 to P7. For each period, the mean water level variation is estimated for each node of the mesh by calculating the difference of the mean water level during this period by the water level at the start of this period.
For the period P4b, with wind directions close to 330°, and high wind speeds (average value 8 m·s −1 ), globally there is a decrease in the water level in the eastern part of the lagoons, and an increase in the western parts. For this period, there is a high transfer of water from the Vaccarès lagoon into the Impériaux and the Lion/Dame lagoons (see Figure 10). For the 21 days simulated, the minimum volume in the Vaccarès lagoon, and the maximum volume in the Lion/Dame lagoon are observed for the period P4b.
At the end of period P4a, the highest water levels are observed in the eastern parts of the lagoons, and the smallest in the western parts (not shown in this paper). This explains why for period P4b (with a transfer of water from the Vaccarès lagoon into the Impériaux and the Lion/Dame lagoons) the highest variations in the Impériaux and Lion/Dame lagoons are observed in the western parts. This highlights the need to have a good knowledge of the initial water levels in the system to assess their variations with wind.
For period P5, there is globally a decrease in the water levels in the Impériaux and the Lion/Dame lagoons, and an increase in the Vaccarès lagoon.
For period P6, there is a transfer from the Lion/Dame lagoon into the Vaccarès lagoon, and a transfer of the Vaccarès lagoon into the Impériaux lagoon. In each lagoon, there are negative variations in the southeastern part of the lagoons, and positive variations in the northwestern parts. For the 21 days simulated, the minimum volume in the Lion/Dame lagoon is observed for the period P6, see Figure 10d.
For period P7, there is mainly a transfer of water from the Impériaux lagoon into the Vaccarès lagoon, and from the Vaccarès lagoon into the Lion/Dame lagoon, explaining the positive mean variations observed in the Lion/Dame lagoon and in the eastern parts of the Vaccarès lagoon.
An interesting point concerns the two maximum volumes in the Impériaux lagoon, which are observed for periods P4b and P6. Interestingly, these two maxima are observed for northwesterly winds in one case (period P4b) and for southeasterly winds in the other (period P6). For the northwesterly winds, water flows into the channels C2 and C3 (see Figure 3) are both from the Vaccarès lagoon into the Impériaux lagoon (see Figure 10f, zoom 1), explaining the relatively fast increase in the volume of the Impériaux lagoon. For the southeasterly winds, water flows in C2 are always from the Vaccarès lagoon into the Impériaux lagoon, whereas water flows in C3 are sometimes Mean water levels variations (m) from the Impériaux into the Vaccarès lagoon, and sometimes from the Vaccarès into the Impériaux lagoon (see Figure 10f, zoom 2). This explains the slower increase in the volume of the Impériaux lagoon, compared to period P4b (see Figure 10c).

Discussion
The Vaccarès lagoon system is an interesting shallow coastal lagoon to study wind-driven circulation patterns. It is indeed unusually shallow for its size and is made up of three sub-lagoons with very different geometries, and with several islands.
Our results show that large amounts of water can be exchanged between the Vaccarès and the Impériaux lagoons, compared to the amounts exchanged between the Vaccarès and the Lion/Dame lagoons. This result is expected to have an impact on the biogeochemical dynamics of the entire Vaccarès lagoon system. This impact will be studied in the future.
In the Vaccarès lagoon, the deepest of the three lagoons, circulation patterns generated by a steady wind are characterized by depth-averaged flow velocities in the same direction as the wind along the shallows, and by slower counter-flows that oppose the wind-direction in deeper parts of this lagoon. These different flows are interlinked by circulatory gyres. These observations are in agreement with the works of [31][32][33][34]. These basic features of the circulation patterns have been observed at many other systems worldwide, as can be seen in the Cabras lagoon in Sardinia, Italy [16], in the Mar Menor lagoon in Spain [35], in the Nador lagoon in Morocco [36][37][38], in the southern part of the Curonian lagoon [39], and in the St. Lucia estuarine lake in South Africa [17]. However, the originality of the Vaccarès lagoon system consists in the Impériaux and Lion/Dame lagoons, which are quite unusual due to their shallowness and complex geometries with several islands. This leads to a large number of circulatory gyres (see Figure 8a,b) whose location and scale are not easy to forecast without hydrodynamic simulations. Wind-induced currents are quite unusual in these two lagoons, in comparison with those observed in lagoons with less complex geometries (Vaccarès, Cabras, Mar Menor, Nador lagoons).
Regarding wind-driven circulations, our results illustrate the differences that exist between the different sub-lagoons of the Vaccarès lagoon system, and the differences in the water exchanges between these sub-lagoons. This is expected to have an influence on the mixing processes in this system. As an example, [40] developed a coupled hydrological and chemical model to simulate time-dependent concentrations of the two herbicides oxadiazon and MCPA, and the two subsequent degradation products of MCPA in the three lagoons of the Vaccarès lagoon system. This model used the assumption of complete mixing. It gave good results for the concentrations in the Vaccarès lagoons, but clearly underestimated the concentrations in the Impériaux and the Lion/Dame lagoons. The authors of this study did not succeed in explaining this underestimation. The results given here (see Figures 8a,b  and 9) show that the assumption of complete mixing in each lagoon cannot be made for short time scales, which certainly explains the bad results in the two lagoons not directly in connection with the agricultural drainage channels (FUM and ROQ, see Figure 1).

Conclusions
This work on the Vaccarès lagoons system illustrates how wind can induce complex currents in a large, shallow, multi-basin lagoon system subjected to a high natural hydrological variability due to the combination of a Mediterranean climate and the artificial hydrological regime imposed by human water management. This is the first study to explore the spatial structure of wind-driven circulation patterns for the whole Vaccarès lagoon system.
Several modeling simplifications were used in this study. First, the model considers a spatially constant Strickler coefficient. It would be interesting to acquire additional experimental data in order to conduct a more in-depth study of the spatial variation of this coefficient (characteristics of the bottom, sea grass distribution).
The dependence of water density on temperature and salinity was not taken into account in the simulations. To correct this, future simulations will incorporate salt and temperature dynamics. Another point to test would be the influence of wind-generated surface waves on the currents simulated.
All system boundaries were specified as closed for the simulations, which are shown in Figures 8a,b and 9. The relative magnitudes of forcing by drainage freshwater and sea water inputs, in comparison with wind forcing, have to be studied.
To go further in the assessment of the influence of wind-induced currents on the biogeochemical functioning, it is planned to study the water transport time scales to characterize the hydrodynamic of the three sub-lagoons of the Vaccarès lagoon system (Vaccarès, Impériaux, Lion/Dame), taking into account forcing by freshwater drainage and sea water inputs. See [41] for more details.
This model is the first necessary step with a view to obtaining a greater understanding of the biogeochemical operation of the Vaccarès lagoon system.