A new parameterisation of salinity advection to prevent stratification from running away in a simple estuarine model

1 Université catholique de Louvain, Unité de Génie Civil et Environnemental, Louvain-la-Neuve, Belgium 2 Université catholique de Louvain, Centre for Systems Engineering and Applied Mechanics, Louvain-la-Neuve, Belgium Received: 21 April 2008 – Accepted: 23 April 2008 – Published: 2 June 2008 Correspondence to: S. Blaise (sebastien.blaise@uclouvain.be) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
Estuaries and their regions of freshwater influence (ROFIs) have been studied for a long time. They exhibit strong gradients of several variables: salinity, temperature, plankton and nutrient concentrations can vary over a wide range of values, strongly impacting physical and biological processes. For instance, complex dynamics, influenced 15 by tides and input of freshwater from rivers, have a strong influence on the growth of phytoplankton (Lucas et al., 1998(Lucas et al., , 1999. This work focuses on estuarine dynamics, especially on the evolution of stratification. The latter is a key player in vertical mixing, which influences directly the vertical fluxes of heat, salt, momentum and nutrients (Simpson et al., 1990). Many studies mechanisms were described in detail by Simpson et al. (1990). Several models were applied to simulate and understand the evolution of stratification in estuaries. Linear prescriptive models were first used (Simpson et al., 1991;Nunes Vaz and Simpson, 1994;Scott, 2004). Then, several authors turned to one-dimensional non linear models (Monismith et al., 1996;Monismith and Fong, 1996;Nunes Vaz and Simpson, 1994;5 Lucas et al., 1998). Recently, three-dimensional models were used to simulate estuarine flows (Hetland and Geyer, 2004;Warner et al., 2005).
One-dimensional non linear models can be very useful to understand and predict the evolution of stratification in an estuary. They are light and simple to build. They require a minimal amount of data and parameters. Furthermore, they generate simple 10 results, which permits to easily understand the key processes and quickly establish diagnoses. However, one common failure of these models is the generation of runaway stratification : when the tidal amplitude is low, stratification tends to grow without bound due to an inadequate parameterisation of horizontal salinity advection (Nunes Vaz and Simpson, 1994;Monismith et al., 1996;Warner et al., 2005). This paper shows that 15 simple analytical developments can lead to a new version of the model which keeps stratification under control. It is also seen that, in the long run, the model is insensitive to an unrealistic initial stratification.
Herein we use a one-dimensional finite-element water column model. Such finiteelement models and their advantages were described by Hanert et al. (2006Hanert et al. ( , 2007. 20 As mixing is a key player in the evolution of the stratification (Nunes Vaz and , we use the Mellor and Yamada level 2.5 turbulence closure Yamada, 1974, 1982;Galperin et al., 1988) which is well suited for the prediction of stratification in estuaries (Nunes Vaz and Simpson, 1994). This turbulence closure was recently implemented using the finite-element method for one-dimensional  25 and three-dimensional  models.
The physical setting is described in Sect. 2. Then, in Sect. 3, the model is presented. Two parameterisations of horizontal salinity advection, the classical one and a new one, are introduced in Sect. 4 and it is seen rigorously that the new approach prevents strat-

189
We will study the stratification in an estuary, which is generated by the front between freshwater and salty seawater. This front is of a crucial importance for the dynamics of the estuary, notably for the vertical density gradient. Therefore, we will consider a water column located at C in Fig. 1, at a distance L to the sea limit. We assume that the salinity at a distance L upstream of C is of the order of S r and that the salinity at the 10 sea limit is of the order of S s . We also assume that S r and S s are constants satisfying the following condition: In such a configuration, the water velocity is mainly caused by two processes (Simpson et al., 1990):

15
-The presence of freshwater coming from the river creates a density front with the salty seawater (Fig. 2a). This front induces a circulation, with light freshwater going towards the sea at the surface, and dense water going towards the river near the bottom. Due to the bottom friction, this circulation is reduced near the sea bottom.

20
-The tidal circulation, influenced by the shear stress due to the bottom friction, generates a parabolic-like velocity profile (Fig. 2b). This profile induces a transport of freshwater varying over the water column, leading to stratification. The succession of ebbs/floods generates a Strain-Induced Periodic Stratification (SIPS) regime, which can be described as follows: during falling tide, a stable stratification develops, which is reduced by mixing at the end of the falling tide. During rising tide, the salinity profile is unstable and is quickly mixed over the vertical, leading to a non stratified water column.
The combination of these processes can generate different flow regimes. If the tides 5 are dominant, the SIPS regime prevails. When the effect of the horizontal density gradient becomes important compared to the tidal effect, the tidal mixing is not sufficient to annihilate the stratification; this stratification strengthens during each tidal cycle, inducing a persistent stratification regime (Lucas et al., 1998). The presence of different non-synchronous tidal components, by generating an alternation of spring/neap tides, can lead to a succession of SIPS and permanent stratification periods (Simpson et al., 1990;Sharples and Simpson, 1993;Nunes Vaz and Simpson, 1994;Monismith et al., 1996).

Model description
The model used herein is based on that of Lucas et al. (1998) and Monismith et al. 15 (1996). For the flow under study, the impact on density of temperature variations is negligible compared with those of salinity. Therefore, density is assumed to be a function of salinity only. As in Lucas et al. (1998) and Monismith et al. (1996), a linear equation of state is adopted: 20 where ρ and S are the density and the salinity, whose reference values are denoted ρ 0 and S 0 , respectively; β is the salinity expansion coefficient, which is assumed to be constant.
191 If x is the horizontal coordinate increasing toward the sea, the along-estuary horizontal velocity u(t, z) at location C obeys the following momentum equation: where g, η, z and H are the gravitational acceleration, the sea surface elevation, the vertical coordinate pointing upwards with its origin at the sea surface and the water col-5 umn depth, respectively. A slip condition is imposed at the bottom: the bottom stress is parameterised as a quadratic function of the velocity (e.g. Blumberg and Mellor (1987)). The surface stress is set to zero. The turbulent viscosity ν is calculated by means of the Mellor and Yamada level 2.5 turbulence closure Yamada, 1974, 1982) implemented in its quasi-equilibrium version (Galperin et al., 1988;Deleersnijder and 10 Luyten, 1994). The surface slope due to the barotropic tides can be represented as in which T i is the tide period and U i ,max the maximum velocity for the i -th tidal component. The baroclinic pressure gradient can be divided into two contributions (Lucas et al., 1998;Monismith et al., 1996): a term derived from the horizontal salinity gra-15 dient, gβτz, and a term derived from the surface slope generated by the baroclinic flow, −gβτγ H 2 , where τ is the longitudinal salinity gradient. The latter is assumed to be constant and is estimated as τ=(S s −S r )/(2L). The dimensionless coefficient γ is to be tuned in such a way that the residual transport is zero, i.e. the average over a tidal cycle of the depth-integrated velocity vanishes. 20 The salinity S obeys the equation where the eddy diffusivity λ is obtained from the same turbulence closure model as the eddy viscosity. The surface and bottom salinity fluxes are prescribed to be zero:

Parameterisation of the horizontal salinity advection term
In the equation governing the evolution of the salinity, most authors (Nunes Vaz and 5 Simpson, 1994;Lucas et al., 1998;Monismith et al., 1996;Monismith and Fong, 1996) assumed the horizontal salinity gradient to be a constant that was evaluated as follows: However, this parameterisation has been identified as the cause of the so-called "runaway stratification", a numerical artefact in which stratification increases indefinitely 10 (Warner et al., 2005). The salinity reaches values that are no longer comprised in the interval [S r , S s ], which is unacceptable. By annihilating vertical mixing, this overestimated stratification corrupts the computation of the evolution of velocity and water properties. This numerical artefact is related to the variation of the forcing terms over the water 15 column in the momentum Eq. (3). Figure 3 shows that, when averaged over a tidal cycle, the sum of each forcing term present in (3) decreases linearly with depth. This variation over the vertical will lead to a seaward tidally-averaged velocity greater in the upper part of the water column than near the bottom. With such a distibution of the velocity, it is obvious that the use of a constant salinity gradient in (5) will inevitably 20 lead to a constantly increasing stratification if the counterbalancing diffusion term is not strong enough. Indeed, at falling tide, the advection of freshwater will decrease with depth whereas, at rising tide, the advection of seawater will increase with depth, causing stratification to grow indefinitely. The apparition of the "runaway stratification" artefact can be avoided by using an alternative parameterisation of the horizontal salinity advection term, inspired by the first-order upwind difference scheme: By introducing u + and u − the positive and negative parts of the longitudinal velocity, and by using relation (8), we can rewrite Eq. (5) as If the velocity is directed toward the sea (u>0), the first term in the right-hand side of (10) relaxes the salinity to its river value S r , the relaxation timescale being L/u + . On 10 the other hand, when the velocity is directed toward the river, the salinity is relaxed toward S s with a relaxation timescale equal to L/|u − |.
It is interesting to notice that resorting to this new parameterisation is equivalent to add to the classical formulation (7) a horizontal diffusion term. Indeed, with the parameterization suggested herein, the horizontal salinity advection may be rewritten Clearly, the last term in the equation above may be viewed as the discrete form of the harmonic diffusion operator, the associated diffusivity being |u| L/2. The interpretation of the role of the first two terms in the right-hand side of salinity 20 Eq. (10) suggests that, whatever the horizontal velocity, the salinity should tend to be 194 comprised in the interval [S r , S s ]. In fact, this can be demonstrated rigorously. For an arbitrary large value of t (t→∞), the salinity must obey the following inequalities: implying that stratification cannot grow out of control. We first define the overshooting of the salinity by So, the overshooting is a positive variable that is equal to S(t, z)−S s if the salinity is greater than its sea value S s , and is equal to zero otherwise. Multiplying Eq. (10) by the overshooting and integrating over the height of the water column yields: The manipulations leading to this equation are not trivial, but they are of the same type as those of Appendix C of Deleersnijder et al. (2001). All of the terms in the righthand side of (14) are negative unless the overshooting is zero at every point of the 15 water column. Thus, the quadratic measure of the overshooting tends to zero as time increases, implying that lim t→∞ δ + = 0.
Combining relations (13) and (15) leads to A similar analysis can be performed for the undershooting δ − = max [0, S r −S(t, z)], eventually leading to S(t, z)≥S r for t→∞. Hence, (12) holds valid. QED. Needless to say, it cannot be seen that, when the classical parameterisation (7) is used, the salinity asymptotically remains within the interval [S r , S s ].

5
To illustrate the advantages of the parameterisation designed above, we will simulate the situations described in Sect. 2. All of the simulations are achieved using a timestep of 60 s. The one-dimensional vertical mesh contains 30 nodes. The main physical parameters are similar to those of Nunes Vaz and . The water column depth is 15 m, and the values of S r and S s are respectively 0 psu and 35 psu.
We first consider a SIPS regime similar to that of Nunes Vaz and . There is only one tidal component with a magnitude of U 0,max =1 m/s and a period of T 0 =12 h. The longitudinal constant salinity gradient τ is set to 0.25 psu/km. Figure 4 shows that the SIPS regime is quickly established, with an alternation of stratified/unstratified phases. The tidal mixing at the end of the falling tide is sufficient to an-15 nihilate stratification. The latter is very similar using both parameterisations of salinity advection. However, the constant parameterisation (7) leads to higher peaks of stratification while the latter is limited using the new parameterisation (8). These smaller peaks can be explained by the horizontal diffusion added to the model (11) when we use the new parameterisation of salinity advection. The mean velocity remains rather 20 insensitive to the used parameterisation.
For the SIPS regime simulated above, the two expressions of the salinity advection term led to rather similar results. This is not always the case, especially if a permanently stratified regime is considered, such as that investigated by Nunes Vaz and . Accordingly, the tidal amplitude is decreased (U 0,max =0.5 m/s) to reduce mixing and the longitudinal salinity gradient is increased (τ=0.3 psu/km). All the other parameters remain unchanged. Model results are displayed on Fig. 5. Using the classical 196 parameterisation of the horizontal salinity advection term, the stratification grows out of control to unrealistic values exceeding the imposed bounds, which is the deficiency known as "runaway stratification". As demonstrated in Sect. 4, the stratification remains within the imposed limits when we use the new parameterisation. The slight oscillations show that, even when the stratification is high, it is still influenced by the tides. While 5 the classical parameterisation (7) gives useless results, the new parameterisation (8) gives qualitatively realistic results for a large number of tidal cycles. The spring/neap cycles are now simulated by taking into account two tidal components. The first one has an amplitude of U 0,max =0.8 m/s and a period of T 0 =12.42 h; while the second component has an amplitude of U 1,max =0.46×U 0,max and a period of 10 T 1 =12 h (Nunes Vaz and . Using this combination of tidal components, we generate an alternation of spring and neap tides (Fig. 6). We set the longitudinal constant salinity gradient to the value of τ=0.25 psu/km. It is shown on Fig. 6 that both parameterisations represent a spring-neap cycle of stratification. During neap tides, the stratification grows until the tidal amplitude increases at spring tides. Then, 15 the stratification weakens and comes back to a SIPS regime. However, the classical parameterisation leads to unrealistic peaks of stratification, with salinity exceeding the limits imposed by the river and sea salinities. This is a common phenomenon when using expression (7) (i.e. Nunes Vaz and  in which the difference between bottom and surface density grows during neaps as far as 180 kg/m 3 ). This 20 problem does not occur when the new parameterisation is resorted to.
In this last experiment, we simulate a spring-neap cycles regime giving rise to the runaway stratification artefact. To this aim, the tidal amplitude is decreased to U 0,max =0.7 m/s and U 1,max =0.46×U 0,max , while the longitudinal constant salinity gradient is increased to τ=0.3 psu/km. Figure 7 shows that the classical parameterisation of 25 the horizontal salinity advection term leads to a stratification which increases unboundedly and then cannot come back to the SIPS regime. During successive tidal cycles, the stratification strengthen to excessively large values. The new parameterisation, by limiting the peak of stratification to acceptable values, permits to come back to the 197 SIPS regime during spring tides which is believed to be consistent with observation (Simpson et al., 1990;Sharples and Simpson, 1993).

Discussion
We now investigate the impact of the initial conditions, in particular the stratification prescribed at the initial instant. Figure 8a shows the evolution of the stratification using 5 different initial stratifications for the SIPS regime. For each parameterisation, the mixing is able to annihilate the stratification, yielding a SIPS regime. However, the decrease of the stratification is much faster using the new parameterisation. It was demonstrated in Sect. 4 that, even if we have an overshooting or an undershooting in the initial salinity, this excess will be eliminated by the new parameterisation of the horizontal salinity advection. If the stratification exceeds the upper limit, it cannot strengthen anymore when the new parameterisation is used, whereas the classical parameterisation still generates cycles of increase/decrease of stratification. In a persistent stratification regime (Fig. 8b), using the new parameterisation, the stratification decreases under its upper limit value, and then reaches a regime solution. The solution converges for 15 any initial stratification. This confirms that any overshooting is directly eliminated by that parameterisation. The classical parameterisation, on the other hand, generates a runaway stratification.
Although the new parameterisation of the horizontal salinity gradient has only been applied to the advection term of the salinity equation, it is interesting to know if introduc-20 ing this parameterisation in the momentum equation has a significant impact. Figure 9 shows the evolution of the stratification for a spring/neap cycle using the two different parameterisations of the horizontal salinity for the momentum equations. To avoid the development of runaway stratification, the new parameterisation is used in the salinity equation. We can see that the influence of the parameterisation of horizontal salinity 25 in the momentum equation has a small influence on the stratification, whereas this parameterisation has a strong influence when applied to the advection term of the salinity 198 equation (Fig. 6).
By slightly modifying the equations, the present model could also be applied to the simulation of the tidal straining in a Region of Freshwater Influence (ROFI), for which the stratification induced by a gradient of density is also a key process (Visser et al., 1994). Two components of velocity u and v would be used (Simpson et al., 2002). 5 The new parameterisation of salinity advection would be able to avoid the generation of runaway stratification in a ROFI model, for which this numerical artefact can also occur.

Conclusions
Using simple mathematical developments, a new expression of the horizontal salinity 10 advection term was developed in order to avoid the numerical artefact known as runaway stratification. This method allows for the simulation of rather realistic flows such as spring/neap cycles without any unrealistic stratification peak. It is guaranteed that no over-or under-shooting will be generated and that any initial over-or under-shooting will progressively disappear. The mathematical method we had recourse to for establish- 15 ing the properties of the new parameterisation of horizontal salinity advection may be applied to a wide range of partial differential problems in order to derive a priori upper or lower bounds of their solution. This technique is inspired by Lewandowski (1997). To the best of our knowledge, it has been used in a small number of oceanographic studies only (Deleersnijder et al., 2001;Legrand et al., 2006;Gourgue et al., 2007). In 20 our opinion, that this approach is so rarely resorted to in the realm of oceanography is somewhat regrettable. The present study shows how simple and powerful this method is. funded by the Communauté Francaise de Belgique, as Actions de Recherche Concertées, under contract ARC 04/09-316 and the project "Tracing and Integrated Modelling of Natural and Anthropogenic Effects on Hydrosystems" (TIMOTHY), an Interuniversity Attraction Pole (IAP6.13) funded by the Belgian Federal Science Policy Office (BELSPO). This work is a contribution to the development of SLIM, the Second-generation Louvain-la-Neuve Ice-ocean Model Fig. 1. Physical setting: the stratification is to be simulated in the water column located at point C. The latter is located in a region of high salinity gradient. Its order of magnitude is (Ss − Sr)/2L, where Ss and Sr denote the downstream and the upstream salinity, respectively.

Fig. 1.
Physical setting: the stratification is to be simulated in the water column located at point C. The latter is located in a region of high salinity gradient. Its order of magnitude is (S s −S r )/2L, where S s and S r denote the downstream and the upstream salinity, respectively.  (7) new parameterisation (8) Fig. 6. Simulation of the circulation induced by a succession of spring/neap tides: results obtained using the old (7) (dashed curves) and the new (8)    Tidal cycles old parameterisation (7) new parameterisation (8) Fig. 9. Impact of the parameterisation of the horizontal salinity gradient in the momentum equation. Evolution of the stratification (difference between bottom salinity and surface salinity) for the old (7) (dashed curve) and the new (8) (solid curve) parameterisations of the horizontal salinity advection.