Numerical investigation of algebraic oceanic turbulent mixing-layer models

In this paper we investigate the finite-time and asymptotic behaviour of algebraic turbulent mixing-layer models by numerical simulation. We compare the performances given by three different settings of the eddy viscosity. We consider Richardson number-based vertical eddy viscosity models. Two of these are classical algebraic turbulence models usually used in numerical simulations of global oceanic circulation, i.e. the Pacanowski–Philander and the Gent models, while the other one is a more recent model (Bennis et al. , 2010) proposed to prevent numerical instabilities generated by physically unstable configurations. The numerical schemes are based on the standard finite element method. We perform some numerical tests for relatively large deviations of realistic initial conditions provided by the Tropical Atmosphere Ocean (TAO) array. These initial conditions correspond to states close to mixing-layer profiles, measured on the Equatorial Pacific region called the West-Pacific Warm Pool. We conclude that mixing-layer profiles could be considered as kinds of “absorbing configurations” in finite time that asymptotically evolve to steady states under the application of negative surface energy fluxes.


Introduction
The mixing layer is located immediately below the ocean surface, and its formation is due to atmospheric-oceanic factors of exchange driven by the stress induced by the winds that generates a strong turbulent mixing dominated by vertical fluxes.The dynamics of mixing layers play an important role in the global oceanic circulation and global climate changes.Indeed, the depth of the mixed layer, which is the upper homogeneous part of the mixing layer with almost constant density, is very important for determining the sea surface temperature (SST) range in oceanic and coastal areas.In addition, the heat stored within the oceanic mixed layer provides a source of heat that drives global variability such as El Niño.The mixed layer also has a deep impact on the evolution of polar ice (Goosse et al., 1999), and it is closely related to different aspects of oceanic bio-systems too.The bottom of the mixed layer corresponds to the top of the pycnocline, a zone of high gradients of density (see Vialard and Delecluse, 1998;Defant, 1936;Lewandowski, 1997 for a physical description of the structure of mixing layers).
The Oceanic General Circulation Models (OGCM) include mixing-layer parameterizations in order to better take into account the influence of the atmosphere-ocean surface interactions (Burchard et al., 2005;Wang et al., 2008).Indeed, they incorporate specific turbulence models for mixing layers, and within these the algebraic ones, which parametrize the turbulent viscosity and diffusion by means of algebraic expressions in terms of the gradient Richardson number.The Richardson number represents the balance between stabilizing buoyancy forces and destabilizing shearing forces.These kinds of models were introduced in the 1980s (Pacanowski and Philander, 1981), and they apply to stratified shear flows that are assumed to have reached a vertical equilibrium after the vertical mixing generated by the wind stress has been restabilized by buoyancy forces.The model proposed by Pacanowski and Philander (1981) was modified in several ways in order to obtain a better fitting with experimental data (Gent, 1991).Another kind of improvement was based upon the parameterization of the vertical profile of turbulent kinetic energy (KPP models; Large et al., 1994).In all these models, only vertical eddy diffusion effects are included.More complex and sophisticated parameterizations Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
of the vertical turbulent mixing, such as the k−ε model, have also been developed (Mellor and Yamada, 1982;Gaspar et al., 1990) and used in the context of OGCM.
The mathematical and numerical analysis of these models technically faces hard difficulties, and to the knowledge of the authors, there are no references addressing this analysis in the last cited models at the present date.In this paper we will focus on the classical algebraic turbulent mixing-layer models of Pacanowski and Philander (1981) and Gent (1991), and on a more recent algebraic turbulent mixing-layer model proposed in Bennis et al. (2010).Indeed, for these three models there exists a mathematical study performed by the authors of the present paper, especially focused on the non-linear asymptotic stability analysis under small perturbations, that applies to negative heat surface fluxes (Chacón-Rebollo et al., 2013a, b).In this paper we extend this study with the numerical investigation of the behaviour of the aforementioned models for relatively large deviations from initial conditions with respect to mixing-layer profiles.The applied initial perturbations admit meaningful physical interpretations: on the one hand, we take into account strong heating or precipitation phenomena at the surface, while on the other hand we consider the cooling or evaporation of surface water.These cases have not been analysed in the above-referenced papers, which deal with initial conditions close to mixing-layer profiles.In addition to studying the behaviour of these models for characteristic times of formation of well-developed surface turbulent mixing layers, we also perform an investigation of their asymptotic behaviour.We conclude that mixinglayer configurations attract the perturbed initial conditions in time scales of the order of several days, in agreement with the physics of the problem, and asymptotically evolve into the theoretical equilibria in time scales of the order of several months.

Setting of algebraic closure models for oceanic turbulent mixing layers
Let z ∈ [−h, 0] be the vertical spatial variable, where h > 0 denotes the thickness of the studied flow that must contain the mixing layer, and t ≥ 0 be the time variable.The mixing layer is assumed to be strongly dominated by vertical fluxes, so that the velocity and density of the fluid are assumed to be horizontally homogeneous.The flow is turbulent and well mixed, so the vertical velocity vanishes.The model affects the statistical mean horizontal velocity (u, v) and the statistical mean density ρ as functions of the variables z and t.In the ocean, the density is a function of the temperature and salinity through a state law, so it is considered as a thermodynamic variable.The Coriolis force is neglected, which is a valid approximation in tropical seas.The variables u, v, ρ satisfy the one-dimensional Reynolds-averaged equations: where ν 1 = a 1 + ν T 1 , ν 2 = a 2 + ν T 2 respectively are the total viscosity and diffusion.Here, a 1 , a 2 are the laminar viscosity and diffusion, and ν T 1 , ν T 2 are the vertical eddy viscosity and diffusivity coefficients.These are assumed to be functions of the gradient Richardson number R, defined as: where g is the gravity constant and ρ r is a reference density for the sea water.
The eddy coefficients corresponding to the Pacanowski-Philander (PP) model are given by: where: with a 1 = 10 −4 , b 1 = 10 −2 , a 2 = 10 −5 (units: m 2 s −1 ).The Gent model is just a variant of the PP model, designed to better fit experimental data, given by: with a 1 = 10 −4 , b 1 = 10 −1 , a 2 = 10 −5 (units: m 2 s −1 ).In general the PP and Gent models become numerically unstable if the initial conditions are physically unstable configurations, corresponding to R < 0 (Bennis et al., 2008).In Bennis et al. (2010), a model of the eddy diffusion that remains numerically stable for a large range of negative gradient Richardson numbers was introduced: with the same constants as the PP model.
In formulas (3) to (5), the eddy coefficients ν T 1 and ν T 2 are defined as functions of the gradient Richardson number through the terms (1 + γ R) n appearing at the denominator.Hereafter, these three formulations will be denoted respectively by R213, R23 and R224, where the integer values are the exponents of (1 + γ R) in the definitions of ν T 1 and ν T 2 .The eddy coefficients defined by relations (3) to (5) all present a singularity for a negative value of the gradient Richardson number.In Fig. 1, we have plotted the curves ν 1 = f 1 (R) and ν 2 = f 2 (R) obtained with the different formulations.In Eqs. ( 3) and ( 4), the diffusivity coefficient ν 2 becomes negative for negative values of R, and therefore these models are no longer valid.
We shall consider the following initial and boundary conditions for problem (1): The Neumann boundary conditions at z = 0 represent the fluxes at the sea surface that model the forcing by the atmosphere.In particular: Q u , Q v are the surface momentum fluxes and Q ρ represents thermodynamic fluxes.The momentum fluxes are given by , where ρ a is the air density, and V u , V v respectively are the stresses exerted by the zonal and the meridional winds: with U a = (u a , v a ) the air velocity at the atmospheric boundary layer, and C D (= 1.2×10 −3 ) a friction coefficient (Kowalik and Murty, 1993).
Note that models (1)-( 6) are not expected to describe all the interaction phenomena occurring in the mixing layer.Its purpose is mainly to give a better understanding of algebraic closure models for oceanic turbulent mixing layers.Therefore, we shall use simplified equations governing the variables u (zonal velocity), v (meridional velocity) and ρ (density).In practice, mixing-layer models are coupled with 3-D models of global oceanic circulation that yield the boundary values for velocity and density at the bottom of the layer.This allows one to use finer grids (in vertical direction) to compute the heat exchange through the ocean surface (Wang et al., 2008).The coupling of 1-D mixing-layer models to 3-D OGCM for the inner oceanic flows means that the 1-D model may be affected by multidimensional perturbations.In Chacón-Rebollo et al. (2013b), the authors tested whether, in tropical seas, replacing a multidimensional model with a 1-D mixing-layer model provides accurate results.

Equilibria states of a perturbed model
We show in this section that there exist steady solutions to problems (1)-( 6) for perturbed data.These perturbations may correspond to errors in the experimental measurements, roundoff computational errors, errors in the boundary data coming from the approximate solution of the 3-D global model, etc.The steady solutions correspond to an equilibrium between destabilizing wind shear effects and stabilizing cooling surface heat fluxes.
Let us consider the perturbed model: with the perturbed boundary and initial conditions: The existence and uniqueness of smooth equilibria for this problem are stated in Chacón-Rebollo et al. (2013b), and are given by the following: Theorem 1 Assume that for any z ∈ I the implicit algebraic equation: where G(z) is the function defined by: admits at least a solution R e .Then, to each solution R e there exists a unique associated equilibrium solution of problem ( 7)-( 8) in [H 1 (I )] 3 , given by: The existence of solutions of the algebraic equation ( 9) is ensured under the following: and for some λ ∈ (0, 1) and i = u, v, ρ it holds: The assumption Q u > 0, Q v > 0 means that the wind velocity acts as a destabilizing agent for the mixing layer flow, while Q ρ < 0 means that the surface heat flux plays a stabilizing role.We conclude that for all considered models, there exist steady solutions to problems (1)-( 6) for perturbed data given by Eq. ( 10).
It can be proved that under Hypothesis 1 there exists a unique gradient Richardson number R e for Bennis et al.
Remark 1 The equilibria of unperturbed model ( 1) are studied in Bennis et al. (2010).In that case, R e does not depend on z, and, as a consequence, the equilibrium profiles for velocity and density are linear.The equilibria for perturbed models ( 7)-( 8) provided by Theorem 1 converge to those of the unperturbed model as the perturbations in the data vanish.

Numerical tests and results
In this section, we study the application of oceanic turbulent mixing-layer models in the equatorial Pacific region, called the West-Pacific Warm Pool.The West-Pacific Warm Pool is an area located at the Equator between 120 • E and 180 • E, where the temperature is high and almost constant along the year (oscillating between 28 and 30 • C).The precipitations are intense and, as a consequence, the salinity is low.In particular, we are interested in perturbing real initial data taken from the TAO array (McPhaden, 1995).These data correspond to profiles close to mixing-layer configurations that have already been analysed by the authors of the present paper (Chacón-Rebollo et al., 2013b;Rubino, 2011).In these previous studies, we have concluded that by starting from initial conditions close to mixing-layer profiles, we reach a well-developed surface turbulent mixing layer within two days, and mathematical stable equilibria within two months approximately.Here, we are mainly interested in analysing whether the formation of a homogeneous mixing layer first, and then of theoretical equilibria, are reachable even starting from initial conditions far away either from mixing-layer profiles and steady-state solutions, within the respective characteristic times.
To numerically solve systems ( 7)-( 8), we discretize the initial boundary value problem by linear piecewise finite elements.To describe the numerical scheme, assume that the interval I = [−h, 0] is divided into N subintervals of length  z = h/N, with nodes z i = −h + i z, i = 0, . . ., N , and construct the finite element space: To discretize the equation for u, for instance, we implement the semi-implicit method: where: and we consider similar discretizations for v and ρ.This discretization, under certain hypotheses, is stable and verifies a maximum principle (Bennis et al., 2010).In Chacón-Rebollo et al. (2013b), we developed a specific numerical analysis for the discrete equilibria of systems ( 7)-( 8), with the actual viscosities corresponding to the PP, Gent and Bennis et al. models, proving the existence, uniqueness and convergence to continuous equilibria.Also, in Chacón-Rebollo et al. (2013a), we proved the non-linear stability of these equilibria, with slight initial perturbations of the steady solutions.
Here, we present two tests.In Test 1, we apply the numerical schemes ( 11)-( 12) with a large negative deviation from realistic density initial conditions in the absence of convection (i.e.∂ z ρ 0 < 0 for any z ∈ (−h, 0)).The perturbation applied in this case intends to simulate strong initial heating or precipitation processes at the surface.In Test 2, we consider an initialization of the code by a large positive surface deviation from a real density profile in the presence of convection (i.e.∂ z ρ 0 > 0 for some z ∈ (−h, 0)).The perturbation applied in this case simulates strong initial physical processes such as the cooling of the surface water or evaporation.All the numerical experiments are grid size-and time step-independent, in the sense that the results remain practically unchanged as z and t decrease.

The initial and boundary data
We perturb initial data available from the TAO array.The TAO project aims at studying the exchange between tropical oceans and the atmosphere, providing data often used in many numerical simulations.In particular, the velocity data come from the acoustic Doppler current profiler (ADCP) measurements.To obtain the actual profiles, we approximate the data measured at 0 • N, 165 • E by a linear interpolation, and we largely perturb the density at the surface, leaving the velocity profiles unchanged.We consider a negative buoyancy flux acting on the sea surface equal to −10 −6 kg m −2 s −1 in all cases, so that the only turbulence source is the wind stress.This negative surface heat flux is in agreement with Gent (1991).The studied layer depth is 100 m.

Test 1
In this test, we consider a perturbation of the initial conditions corresponding to the TAO data for the 15 June 1991.
The initial profiles are displayed in Fig. 2. The initial zonal velocity presents a westward current at the surface and, below it, an eastward undercurrent whose maximum is located at about −55 m.Deepest, we observe a westward undercurrent.The initial meridional velocity presents a southward current whose maximum is located at about −90 m.The initial density profile presents a strong negative deviation from the original data at the surface (first 20 m), in order to move the density profile away from a (homogeneous) mixing-layer configuration.Physically, we are in the absence of convection, since the initial density profile has a negative derivative along the water column, and by the perturbation we are simulating an important initial increase in heating or precipitation phenomena.At the surface, we impose a zonal wind equal to 8.1 m s −1 (eastward wind) and a meridional wind equal to 2.1 m s −1 (northward wind).
The formation of a well-developed mixing layer is achieved by integrating the various models for time t = 192 h (8 days).The grid spacing is z = 5 m and the time step is t = 60 s.The corresponding numerical results are displayed in Fig. 3.We consider hereafter a standard definition of mixed layers (Thomson and Fine, 2003;Peters et al., 1988) that states that the base of the mixed layer is the depth at which the density changes by 0.01 kg m −3 .The final density profile displays a similar mixed layer for the R213, R23 and R224 models of about 20 m, in agreement with the observations reported by Boyer Montégut et al. ( 2004).Furthermore, the pycnocline simulated by the three models is similar.In the upper oceanic layer, the surface currents for the R213 and R224 models show almost the same behaviour, while the R23 model underestimates these currents.We notice an increase in the zonal and meridional surface currents in comparison with the initial profiles, which is in agreement with the application of a northeasterly wind.The surface current behaviour can be explained by the final viscosity and diffusivity values, displayed in Fig. 4 for all models, where in particular we observe that the R23 model produces the strongest viscosity and diffusivity.
To obtain the steady states of the flow, we integrate the various models for t = 10 000 h (about 14 months) with t =  1 h.In Fig. 5, we can notice the formation of linear steady profiles for velocity and density, conforming to the theoretical expectations (Bennis et al., 2010).Similar numerical results are obtained for the R213 and R224 models, while the R23 model diverges from these results in all the water columns.In Fig. 6, we can observe a monotonic numerical convergence to the steady states for all models, the R23 model being the one that reaches a stable equilibrium first (after about 3000 h, i.e. 4 months approximately).Here, we compute the residual values as: and we refer to a stable equilibrium when r n < 10 −6 .Collecting data at time t = 10 000 h implies the consideration of a subsequent relaxation time until r n ∼ 10 −12 is obtained for all models.This shows the non-linear stability of equilibrium solutions that act as point-wise attractors whenever we apply a negative buoyancy flux at the surface, as in this case.In particular, the steady states attract configurations corresponding to mixing layers, which appear as intermediate transient states reached by the asymptotic evolution of the flow to mathematical equilibria.

Test 2
The code is initialized with a perturbation of the 15 November 1990 data (see Fig. 7).The initial zonal velocity profile displays an eastward current all along the water column, whose maximum is located at the sea surface.The initial meridional velocity displays a southward current whose maximum is located at about 40 m of depth.The initial density profile displays a strong positive deviation from the original data at the surface (first 20 m), in order to move the density profile away from a (non-homogeneous) mixing-layer configuration, and to create a deep (first 50 m) significant surface zone of negative initial gradient Richardson numbers.We are thus in the presence of convection, where convection phenomena are in particular linked to physical processes such as the cooling of surface water or evaporation.Note that the presence of static instability zones in the density profiles can be detected with a certain frequency in the TAO data.In Fig. 8, we display the initial gradient Richardson number, which is effectively negative for the first 50 m of depth.This fact is reflected in the initial diffusivities of the various models as shown in Fig. 9, where in particular we observe negative values for models R213 and R23, while for model R224 we always have positive values.Physically, negative diffusivity cannot exist, so we cannot use R213 and R23 models in their original formulations here (see Eqs. 3 and 4), described in Pacanowski and Philander (1981) and Gent (1991).In practice, ocean modellers bypass this problem by extending the eddy diffusivities to those regions with positive constant values.An example is given by a modified Pacanowski-Philander model, developed in Timmermann and Beckmann (2004), which gives good results in the presence of convection.It consists of adding a constant value of 10 −2 m 2 s −1 to the turbulent diffusion (and viscosity), depending on depth.In particular, the modified expressions depend on the Monin-Obukhov length that characterizes the oceanic surface boundary layer where the wind effects are strong, and it is computed from a diagnostic equation given  by Lemke (1987).Timmermann and Beckmann (2004) show that the modified model gives good results in the Weddell Sea, an Antarctic geographical region where important convective phenomena occur, but it needs the computation of supplementary parameters.The R224 model presents the advantage that it can be used without needing to compute any additional parameter, which is an important feature in situations where we have to simulate different turbulent regimes.We thus only use the R224 model for this test.The large value displayed by the initial R224 diffusivity corresponds to a negative gradient Richardson number near its singularity, that is at R = −0.2.
The formation of a homogeneous mixing layer is achieved by integrating the R224 model for t = 48 h (2 days).The grid spacing is z = 5 m and the time step is t = 60 s.We impose a zonal wind equal to 7.3 m s −1 (eastward wind) and a meridional wind equal to 2.1 m s −1 (northward wind) at the surface.The results are displayed in Fig. 10.We observe that the application of a negative buoyancy flux at the surface permits us to stabilize the water column, producing a 70 m deep homogeneous mixed layer.We notice an increase in the zonal and meridional surface currents, which is in agreement with the application of a northeasterly wind.During the computation, the model is able to automatically adjust the diffusivity, reducing the initial peak until converging to a range between 10 −4 m 2 s −1 and 10 −2 m 2 s −1 , which is in agreement with Osborn and Cox (1972).By estimating the diffusivity with measurements of very small scale vertical structures, they have shown that it fits in the range [10 −6 m 2 s −1 ; 10 −1 m 2 s −1 ] in the studied region, i.e. the West Pacific Warm Pool.
To obtain the steady states of the flow, we propose to integrate model R224 for t = 10 000 h (about 14 months), with t = 1 h.As the equilibrium solutions are unique in the case of model R224, we can compute the theoretical solutions and compare them with the numerical ones.In Fig. 11, we observe that we obtain practically identical linear profiles for the theoretical and the numerical solutions.In Fig. 12, the residual values display a good monotonic numerical convergence to the steady states.Thus, although physically stable equilibria correspond to positive gradient Richardson numbers, even initial conditions that present strong and deep vertical instabilities could asymptotically converge to those equilibria, i.e. these initial conditions too are in the range of attraction.

Conclusions
In the numerical experiments, we have observed the formation of mixing layers as well as the reaching of equilibrium profiles by starting from physically reasonable initial conditions far away from these two different configurations, remaining in the respective characteristic times of formation.We have noticed an increase in the zonal surface current by applying an eastward wind at the surface (u a > 0).At the same time, an efficient northward wind at the surface (v a > 0) causes an increase in the meridional surface current.These results are in agreement with the physical reality.Furthermore, looking at the residual values, we observe a correct convergence to steady states, corresponding to rather high levels of turbulent diffusions due to the dissipative nature of the equations of the models.We have also noticed that the convergence to steady states takes time of the order of several months.This time scale is much larger than that needed to generate a typical homogeneous mixed layer, which is of the order of a few days, by applying steady boundary data.This may explain why these equilibria profiles are not found in real surface layers, as the boundary data usually change in time scales of the order of hours.Even the formation of a homogeneous mixed layer is not possible if these data change too fast.From a numerical point of view, the investigation of the asymptotic behaviour of algebraic turbulence models for oceanic mixing layers has been useful on the one hand for showing the effective excellent stability properties of these schemes for negative heat fluxes, and, on the other, for proving that the equilibria states behave as pointwise attractors for intermediate transient mixing-layer configurations, for which the analysed models provide accurate physical predictions.In their turn, mixing-layer profiles act as "absorbing configurations" with characteristic times of a few days.
where D u , D v and D ρ are functions of z, and σ u , σ v , σ ρ , U b , V b and R b are constants.Let us denote I = (−h, 0), and define the functions

(
2010) model (5) only(Chacón-Rebollo et al., 2013b), which could be interpreted as the intersection of the curves h z