Modeling confinement in Etang de Thau: numerical simulations and multi-scale aspects

The main purpose of this paper is to implement a mathematical model of confinement in a realistic configuration; the model we consider here relies on a previous work by two of the authors in [FR12]. The concept of confinement was introduced by Guelorget and Perthuisot [GP83a] in 1983. It has latter been widely used, studied, discussed and tested (see Guelorget and Perthuisot [GP83b, GP83a], Guelorget, Frisoni and Perthuisot [GFP83], Guelorget et al. [GGLP90], Ibrahim et al. [IGF+85], Debenay, Perthuisot and Colleuil [DPC93], Redois and Debenay [RD96], Barnes [Bar94], Frenod and Goubert [FG07] and Tagliapietra et al. [TSG09]) leading to the conclusion that it is a pertinent parameter controlling the features of living benthic population in paralic ecosystems. Benthic species are species living on the seabed and paralic ecosystems are ecosystems encountered in estuaries, lagoons and closed bays. This is a first motivation to conduct numerical simulations of confinement in a lagoon.


Introduction
The main purpose of this paper is to implement a mathematical model of confinement in a realistic configuration; the model we consider here relies on a previous work by two of the authors in [FR12]. The concept of confinement was introduced by Guelorget and Perthuisot [GP83a] in 1983. It has latter been widely used, studied, discussed and tested (see Guélorget and Perthuisot [GP83b,GP83a], Guélorget, Frisoni and Perthuisot [GFP83], Guélorget et al. [GGLP90], Ibrahim et al. [IGF + 85], Debenay, Perthuisot and Colleuil [DPC93], Redois and Debenay [RD96], Barnes [Bar94], Frénod and Goubert [FG07] and Tagliapietra et al. [TSG09]) leading to the conclusion that it is a pertinent parameter controlling the features of living benthic population in paralic ecosystems. Benthic species are species living on the seabed and paralic ecosystems are ecosystems encountered in estuaries, lagoons and closed bays. This is a first motivation to conduct numerical simulations of confinement in a lagoon. In addition, since the confinement of an area is closely related to the amount of available nutrient, we will implement our model in a lagoon where shellfishes are farmed, here Étang de Thau, located on the French Mediterranean shoreline.
We start with recalling in Section 2 the mathematical model introduced in [FR12]. Then, we perform in Section 3 the numerical simulations of confinement in Étang de Thau, and discuss how to get rid of the spurious oscillations that occur because of the lack of smoothness of the computational domain. Finally, we introduce in Section 4 a method that allows to account for multi-scale aspects of embedded lagoons. This method is based on the theory of absorbing boundary conditions first introduced in [EM77], and is implemented here in a very simple case.

Mathematical modeling
In this section, we want to introduce the mathematical model that will be used in the sequel to compute confinement fields in Étang de Thau. This model is similar to the PDE model introduced by two of the autors and we refer the interested reader to [FR12] for additional details.

Definitions
The first mathematical definition of confinement was provided by Frénod and Goubert (see [FG07]) as a controlling parameter of tide-influenced paralic ecosystems: Definition 1 The confinement value at any point of the lagoon is the time for the sea-water to reach this point.
In order to account for possible time-oscillations (i.e. tidal oscillations), Frénod and Rousseau extended the definition in [FR12]: The instantaneous confinement is, at a given point of the lagoon and at a given time, the amount of time the water which is at the considered time at the considered point has spent inside the lagoon water mass.
Definition 3 The effective confinement is the time-average of the instantaneous confinement over its oscillating period.

Transport equation for confinement
We consider (see Figure 1) that the lagoon is a cylinder with base a regular, connected and bounded domain Ω ⊂ R 2 with boundary ∂Ω. This boundary is shared into Γ and Γ 0 with Γ ∩ Γ 0 = ∅. Any point in Ω is denoted (x, y). The lagoon seabed is described by a piecewise continuous function b : Ω −→ R + , where b(x, y) represents the bathymetry level at the horizontal position (x, y) ∈ Ω. The water altitude h is such that h > sup Ω {b}, exluding outcrops. In summary, the geometrical model of the lagoon writes:  In order to compute the instantaneous confinement, we use a passive tracer g t (see below) advected by the water velocity field u u u. As shown in [FR12], this model is compatible with any lagoon geometry (shape and bathymetry), with the only restriction that intertidal zones and seabed outcrops are not taken into account. The idea developped in [FR12] is to solve the following transport problem: for any time t > 0 and given a sufficiently large time T , the solution is such that g t (T, x, y) is the value of the instantaneous confinement at time t ∈ R + and position (x, y) ∈ Ω, where the water velocity field u u u(t, x, y) may be computed by solving the following equation: where F is a function defined on Γ such that : The velocity field u u u can be separated in several "elementary" velocity fields, each of those being solely induced by one single process (such as evaporation, tide, river input, etc.). Consequently, depending on those processes, the function θ can model several phenomena.
Remark 1 If we consider a tide-submitted lagoon, the domain Ω is time-depending. The reader can refer to [FR12] for more details.

Numerical simulations in a realistic geometry
In order to extend the results of [FR12], we now proceed to a numerical simulation of confinement in a realistic geometry. We consider the case of a lagoon on the French Mediterranean shoreline: Étang de Thau (see Figure 2).

Étang de Thau
Étang de Thau is one of the largest lagoons on the French Mediterranean coastline. It is more than 20km long, from Balaruc-les-Bains (NE end) to Marseillan (SW). A lot of shellfishes live in this lagoon: Thau is famous for its oyster and mussel farms (see [DDG75]). This lagoon is constantly monitored because it is a very sensitive area with regard to eutrophication [AAB + 99]. The sale of oysters and mussels has even been prohibited several times since the late 80's. Indeed, eutrophication induces an ecosystem deficiency: high densities of nitrogen and phosphor make macroalgae and phytoplankton proliferate, which disturbs herbarium development. Finally, this leads to spats infirmity in the shellfish farms and ends with weak crops. In order to quantize eutrophization of a lagoon, we first consider its confinement, with the natural idea that the more an area is confined, the more it is subject to eutrophication.
In Étang de Thau, most of the seawater flux comes from the Mediterranean Sea by the Graus de Pisse-Saumes (0.75 × 10 6 m 3 by day) and Sète's channels (3.7 × 10 6 m 3 by day). Fresh water comes from the rain (48 × 10 6 m 3 by year) and river input (30 × 10 6 m 3 by year) [A.62]. In our numerical simulations, we will account for both seawater sources, in the (simple) case where the evaporation rate θ > 0 is constant, letting time and space variations of θ to subsequent studies. Thanks to data provided by the Languedoc-Roussillon Region, we have built two meshes of the lagoon including its bathymetry map (see Figure 2) on which we will perform numerical simulations of the model introduced in Section 2.

Viscosity and mesh influences
In order to compute confinement in Étang de Thau, we plugged the mesh and bathymetry informations into the software implemented by Frénod and Rousseau in [FR12], thanks to the finite element toolbox FreeFem++ [HPLH04]. Because of non regular boundaries, the corresponding numerical solutions include spurious numerical oscillations. To overcome this problem, we could either add some artificial viscosity in the model, or try to use a finer mesh. In the case where a viscous model is used, we have: where ν is the (dimensionnal) artificial viscosity. Figures 3c and 3d correspond to the same coarse grid and the simulation is run with the same parameters in Equation (5), except the viscosity that is changed from 100 to 400. 1 The reader can see the oscillations occuring in Figures 3a and 3b (white parts of the plot), close to the two entrances; these oscillations disappear with an increased viscosity (Figures 3e and 3f), but the solution has artificially spread inside the domain (Figures 3c and 3d).
Alternatively, we tried to avoid the numerical oscillations whilst keeping a small viscosity: this necessitates a finer mesh, and was done in Figure 4d. The numerical simulation is improved with respect to Figure 3c, without any artificial spread (as in Figure 3d). Naturally, this last simulation is computationally more demanding, since the time step has to be reduced together with the mesh size (for obvious stability reasons).

Multi-scale domain decomposition for embedded lagoons
In this section, we want to account for the possible multi-scale aspect of confinement. Indeed we considered here Étang de Thau as a lagoon of the Mediterranean Sea, but one could also focus on a small part of Thau and consider it as one of its sub-lagoons. On the other side, the Mediterranean Sea can be seen as a (large) lagoon of the North-Atlantic Ocean, etc. There is thus a need to be able to compute a domain confinement while considering this multiscale aspect. In the sequel, we will consider a reference lagoon Ω and perform a classical confinement computation (see above), providing C Ω in Ω. Then, we will truncate the original domain in order to remove a small part ω of it, corresponding to a sub-lagoon, and compute the solution C Ω\ω on the truncated domain Ω \ ω. If the new boundary conditions that are required on Ω \ ω are well chosen, we will see that the two numerical solutions C Ω and C Ω\ω match in Ω \ ω.
The new boundary that is introduced when we remove ω from Ω has to be carefully studied. Indeed, we should implement absorbing boundary conditions (ABC) on it (see [EM77,HS89]). In practice, ABC are very difficult to find and implement; indeed, they lead to nonlocal boundary operators (see [EM77]) that are not numerically suitable. But in our case, the velocity equation is very simple, and the corresponding ABC are no-slip boundary conditions: where u ext n u ext n u ext n is chosen accordingly with the volume of ω (the truncated part of the domain). These no-slip conditions for u u u write as non-homogeneous Neumann boundary conditions for the velocity potential (see [FR12]). As far as the confinement field is concerned, we choose to impose an homogeneous Neumann condition on the new boundary. This is the simplest approximation of the ABC that one can use for advection-diffusion equations (see [Hal86]). We will consider improved boundary conditions (in the sense of [Hal86]), namely first order or second order conditions, in future works. Left: confinement computed in the truncated domain Ω\ω with homogeneous Neumann boundary conditions. Right: relative error between solutions computed in Ω and in Ω \ ω.
In Figure 5a we plot the confinement computed in the truncated domain Ω\ω: it matches almost perfectly the previous simulations (see Figure 4d) in Ω \ ω. Actually, one can see in Figure 5b that there is a slight mismatch located by the new boundary (≈ 10% in L ∞ relative error, 0.6% in L 2 relative error), which is due to the approximation of the ABC (see discussion above).

Conclusion
We achieved two goals in this paper. On the one hand, we confirmed that the mathematical model introduced in [FR12] is suitable for the simulation of confinement in realistic lagoons such as Étang de Thau; the model may now be used and validated with realistic data (i.e. time series for the function θ, see Equation (3)). On the other hand, we introduced in Section 4 a method to perform multi-scale simulations of confinement in the case of embedded lagoons. We will consider theoretical and numerical improvements of this method (boundary conditions, use of homogenization methods) in future works.