Paralic confinement computations in coastal environment with interlocked areas

This paper is in the continuity of a work program. Its goal is to develop an approach of the paralic confinement usable from the modeling slant, before implementing it in numerical tools. More specifically, we here deal with the multiscale aspect of the confinement. If a paralic environment is separated into two (or more) connected areas, we will show that is possible to split the confinement problem into two related problems, one for each area. At the end of this paper, we will focus on the importance of the interface length between the two subdomains.


Introduction
Paralic confinement is one of the most pertinent parameters controlling the features of living species in paralic environments (i.e. environments such as lagoons, estuaries, bays, etc.). It was first introduced by Guélorget & Perthuisot [GP83b,GP83a] and it is linked with nutrient concentration of water in the paralic environments. It was widely discussed and tested in Guélorget [RD96], Barnes [Bar94], Frénod & Goubert [FG07] and Tagliapietra et al. [TSG09]. Its knowledge in a given paralic environment is an important factor for supporting decision of decision makers acting on the paralic environment. For instance, it can be used to help and choose the exact localization of shellfish farms, or to estimate the impact of the building of a new dyke or dam.
Since the recent works of Frénod & Goubert [FG07], Frénod & Rousseau [FR13] and Bernard, Frénod & Rousseau [BFRss], we know that it is possible to develop a methodology to simulate numerically the confinement in any paralic environment. We now enter a phase of our work program which long term objective is to provide an operational tool to compute the paralic confinement in any point of any paralic environment on earth, only from bathymetry and oceanographic data.
Many questions need to be reached before achieving such an objective and we start here by tackling a modest (but important) one, which is related to the capability of computing separately paralic confinement in two connected areas of a given paralic environment. We shall particularly focus on the interface boundary condition that is required for such a coupling.

Interlocked areas
Coastal environment is made of interlocked areas. For instance, if we take a look at a marsh in a Mediterranean lagoon, we face with the following cascade of areas: Atlantic Ocean -Gibraltar Strait -Mediterranean sea -lagoon entrance -the lagoon -marsh entrance -the marsh. Beside this, coastal environments may present a wide range of scales. For instance in the previously evoked cascade, the marsh is several tens of meters large, while the lagoon size is about ten kilometers. Those two scales are small when compared with the characteristic size of the Mediterranean sea, which is itself small with respect to the Atlantic Ocean dimension. Sizes of transition inlets -Gibraltar Strait and the lagoon entrance -need also to be taken into account.
As it will be recalled in the sequel (see also [FR13]), the numerical computation of the paralic confinement in coastal environment first requires the computation of the water flow essentially going from the ocean to the coastal environment far end, induced by the combined effect of evaporation, tide and fresh water inputs from the rivers. Once the flow is known, a tracer following this flow and undergoing diffusion, is then computed. At the end of the process, this tracer provides the value of the paralic confinement. Because of the wide range of scales appearing in coastal environments (see above), a confinement simulation may rely on several mathematical models that one needs to couple. In the case where the two coupled models are identical, we face a domain decomposition problem. Even if the coupling between the two (ore more) subdomains should actually be two-way, we will focus on one-way exchanges (from the open deep sea to the lagoon, and finally to the marsh). In other words, when we decompose a computational domain in two parts, we will consider a main lagoon (the part that is directly connected to the open sea) and a secondary lagoon (see Figure  1 below). The confinement field in the main lagoon will have to be computed accurately in a truncated domain Ω main = Ω \ Ω seg , while the simulation in the secondary lagoon Ω seg is nothing but another "classical" (i.e. monodomain) confinement simulation, with the main lagoon playing the role of the open deep sea. In the next subsection we recall the model previously developed to compute the paralic confinement field in a lagoon. In subsection 3.3, we introduce the domain decomposition Ω = Ω main ∪Ω seg and the related issues. We pay a particular attention to the interface (and the related boundary conditions) between the main and secondary lagoons.

Domain and equations for paralic confinement computation in a lagoon
We consider (see Figure 2) a lagoon that is a cylinder with base a regular, connected and bounded domain Ω ⊂ R 2 with boundary ∂Ω. This boundary is shared into Γ in and Γ 0 with Γ in ∩ Γ 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: Lagoon = (x, y, z), (x, y) ∈ Ω, b(x, y) < z < h . (1) The upper part is denoted Ω seg . Left: large interface (δ = 20r 1 /100), where r 1 is the radius of Ω main . Right: tinier interface (δ = 10r 1 /100).
Simulations will be conducted with δ/r 1 = 5%, 10%, 15% and 20%.  The instantaneous confinement is -at a given time t and a given position (x, y) -the amount of time the water particle located at the position (x, y) at time t has spent inside the lagoon water mass. It is related to nutrient concentration at position (x, y) and time t. It is the result of two phenomena. In the first place, evaporation process generates a flow from the ocean to the lagoon far end. Essentially, ocean waters can be seen as having a very high nutrient concentration and when those waters travel towards the lagoon, they meet living species that take off nutrients. Then along the travel, their nutrient concentration decreases. Secondly, the nutrient concentration undergoes diffusion because of water molecular and eddy viscosities, but more pregnantly, because of the chemical precesses involved in the dissolution phenomena. This diffusivity is small, when compared with the average consequence of the flow. Nevertheless, it is important in places where the velocity of the flow is small, especially in the lagoon far end. Consequently, in order to compute the instantaneous confinement, we use a passive tracer g t advected by the water velocity field u u u and undergoing diffusion. As shown in Frénod & Rousseau [FR13], 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 [BFRss] is to solve the following advection-diffusion problem: for any time t > 0 and given a sufficiently large time T , the solution g t = g t (τ, x, y) of is such that g t (T, x, y) is a good approximation of the value of the instantaneous confinement at time t ∈ R + and position (x, y) ∈ Ω. Here τ is a variable related to the time to spend into the lagoon and ν is the small (when compared with the average value of u u u) diffusivity coefficient that models the nutrient hability to spread out the water. In system (2), the water velocity field u u u(t, x, y) may be induced by several phenomena (such as evaporation, tide, river input, etc.), which are modeled by the generic function θ. We may compute this field by solving the following equation: where n stands for the unitary vector orthogonal to ∂Ω pointing outside Ω and where F in is a function defined on Γ in such that : We remark 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 Consequently, depending on those processes, the function θ can model one or several phenomena. The sytem of equations (3) is solved thanks to its corresponding velocity potential formulation (see [FR13]). Provided that ∇ × u = 0, we write u = ∇ψ and solve the following Laplace equation for ψ: In subsection 3.3 and section 4 below we will consider equations for the velocity in truncated regions of the lagoon (e.g. equations (10) and (6)). Naturally, these equations will always be solved thanks to the potential formulation, even if this is not explicitly specified.

Domain decomposition for a lagoon with a secondary lagoon
When in a lagoon a clearly separated entity -so called secondary lagoon -exists, we want to split problem (2), (3) and (4) into two problems -a first one set in the secondary lagoon and another one set in the remainder of the lagoon -being connected by conditions on the secondary lagoon entrance. Naturally, we want the concatenation of the results to approximate a solution of the system (2), (3) and (4) set in the whole lagoon with a good accuracy. The way to account for this situation, in what concerns the geometrical aspects, consists in sharing lagoon Ω into three parts (see the left picture in Figure 2): the secondary lagoon Ω seg , the main part of the lagoon Ω main and their common boundary Γ. They are such that Ω seg and Ω main are open subsets of Ω, Γ = ∂Ω seg ∩ ∂Ω main and Ω seg ∪ Ω main ∪ Γ = Ω. Moreover we assume that the secondary lagoon is not near the lagoon entrance, which is translated by ∂Ω seg ∩Γ in = ∅.
The key problem -and the most difficult one -is the obtention of the water flow and of the tracer within the main part of the lagoon without computing them in the secondary lagoon. We now focus on this issue and, in a first place, we build the system of equations to compute the velocity of the water flow in the main part of the lagoon Ω main . The two first equations of problem (3) are retained, because they describe the physics of water transport. The third equality of (3) that translates that water cannot escape from the secondary lagoon through the shore will also be kept. Thirdly, a condition on the interface -that we will discuss hereafter -will be written. Hence, in the main part of the lagoon Ω main , we write the following problem that allows us to obtain the velocity where vector n and function F in have the same definitions as in system (3) and equality (4), where n trans stands for the unitary vector, orthogonal to Γ, and pointing inside Ω main (or outside Ω seg ). The function F is to be determined. It is clear that if the solution u u u of (3) were known, we would choose F = u u u · n trans and then we would obtain a solution u u u main of (10) that would be such that u u u main = u u u |Ω main . Yet, we work under the assumption that the solution of (10) is not known (we indeed want to compute u u u main to have the value of u u u in Ω main ). Anyway, using the Laplace-Neumann compatibility condition (the same that brings us to write (4) for system (3) for u u u), we know an information on F which is and that translates that the quantity of water entering Ω seg -and so leaving Ω main -through Γ compensates for what is consumed by the process modeled by θ over Ω seg . Knowing this information, we can consider that the missing information is the profile (or the shape) of F along the interface Γ. Approximations of this profile, that are classical and known as giving proper results, can be used. For instance, we can chose F as being constant along Γ or being a Poiseuille profile (see [Poi44]).
Having u u u main on hand, and then considering that we consequently know u u u on Ω main with a good accuracy, in order to compute the passive tracer given by (2) only on Ω main , we will consider the following problem: which straightforwardly comes from (2) replacing u u u by u u u main . This system has to be coupled with boundary conditions on the interface Γ. As the diffusivity coefficient ν is small, following Halpern [Hal86], choosing the following Neumann condition will give a solution g main t which will correctly approach g t over Ω main .
Once brought a way to tackle the key problem, we can notice that we can implement a way to obtain the water flow and the tracer within the secondary lagoon. As a decision was made concerning the profile of F , and so concerning F on Γ, we can write the following system to obtain the velocity field u u u seg in the secondary lagoon Ω seg : where n has the same definition as in system (3) and where n trans and F are the one set to solve system (6). On the other hand, g main t on interface Γ can be computed as a result of system (8) and gives the value of the tracer on this interface. This function can be used as a Dirichlet boundary condition for the problem giving g seg t , which will be close to g t in Ω seg . This problem reads: ∂g seg t ∂τ (τ, x, y) + u u u seg (t − T + τ, x, y) · ∇g seg t (τ, x, y) − ν∆g seg t (τ, x, y) = 0, with g seg t (τ, x, y) = g main (τ, x, y), ∀0 < τ < T, ∀(x, y) ∈ Γ.

Numerical simulations
In this section, we present numerical simulations of the lagoon described in Figure 1 (see Section 3). These simulations were performed with the finite element method implemented in the FreeFem++ software [HPLH04]. As in [BFRss] the velocity equation is solved thanks to a Laplace equation on the velocity potential (u u u = ∇ψ and we use P 2 elements for ψ), whereas the advection-diffusion equation on g t is solved thanks to P 1 elements. We consider four different configurations, in order to enhance the importance of the interface width |Γ|. The lagoon Ω is equally split in two parts Ω main and Ω seg , which is the most general (unfavorable) case. Indeed, in cases where |Ω seg | < |Ω main |, the truncation error obtained in the numerical simulations is lower. 1 For each of these configurations, we perform a numerical simulation of confinement in the whole domain Ω, thanks to equations (2) and (3). The corresponding numerical solutions will be considered as reference solutions and denoted (u u u ref , g ref t ).

Simulations without interface information
We now consider the numerical simulation of confinement in the truncated domain Ω main (see Figure 1), in which we look for u u u main , g main t solutions of systems (6) and (8). For the sake of simplicity, we consider a flat bottom and set h − b ≡ 1, so that the boundary function F is such that: where f pr denotes the (unknown) profile of the velocity along the interface Γ and is such that x, y)dσ = 1. In the simulations below, we use a Poiseuille profile for f pr .
We enumerate in Table 1 the L ∞ norm of the relative error between the reference solution g ref t restricted to the subdomain Ω main and g main t (computed in Ω main from (8)- (9)). This error is defined by We observe from the results in Table 1 that the longer the interface, the larger the error. This is due to the lack of information we have at the interface (in particular on the velocity profile, see Section 3 and discussion above).

Simulations with interface information
We now reproduce the numerical simulations of Section 4.1, but instead of choosing a Poiseuille profile for the function f pr , we use the exact profile provided by the knowledge of u u u ref on the interface 2 : u u u main · n trans = u u u ref (t, x, y) · n trans ∀t ∈ R, ∀(x, y) ∈ Γ.
Remark 1 We obviously have that Γ u u u ref (t, x, y) · n trans dσ = Ω seg θ(t, x, y) dx dy, which means that the lack of information on Γ in Equation (13) concerns the velocity profile f pr rather than its average amplitude over the interface.
Then, thanks to the well-posedness of the corresponding velocity equations, uniqueness immediately insures that u u u main = u u u ref |Ω main , that is to say the knowledge of the velocity profile along the interface Γ insures the complete knowledge of the velocity in Ω main . As above, we enumerate the corresponding errors in Table 2, and as expected the errors also depend (in the same manner) on the interface width, but are notably lower than those of Table 1  Furthermore, we can illustrate the quality of the Neumann boundary conditions (9) used for the confinement equation with regard to the confinement diffusivity. Table 3 illustrates that when the diffusivity is low, Neumann boundary conditions efficiently approximate the exact transparent boundary conditions, which was already proved in [Hal86].
Diffusivity ν L ∞ relative error 1.10 −1 0.010308 5.10 −2 0.00323477 1.10 −2 0.00267697 5.10 −3 0.00195359 Table 3: Relative error between confinements g main t and g ref t as a function of the diffusivity ν, in configuration 4 (worst case, large interface width δ = 20r 1 /100). The error increases with the diffusivity. The figure only illustrates the results indicated in the table: the plot corresponds to the infinite relative error (column 2) as the function of the diffusivity (column 1).

Conclusion
In this paper we are interested in the truncation of computational domains in confinement models. This is a very important issue both for classical domain decomposition problems and for the numerical simulation of confinement in limited areas of large lagoons, which we consider here. As soon as the domain truncation is done, the most important question to address is the search for artificial boundary conditions at the new boundary. It is known that their nature strongly depends on the PDE model that drives the considered process. Starting from the confinement model introduced in [FR13], we introduced some boundary conditions in order to limit the numerical error induced by the domain truncation. The chosen confinement condition is a classical homogeneous Neumann boundary condition, which is known to be accurate for small diffusivity values (see [Hal86]). The truncation error is actually mainly due to the lack of knowledge of the velocity profile across the artificial boundary, particularly in the case where the interface and/or the secondary lagoon are large. It would be very interesting to evaluate how some partial informations on this velocity profile (provided by measurements) would improve the corresponding truncation error. We leave this to subsequent studies.