A comparison of three turbulence models with an application to the West Pacific Warm Pool

arXiv:math-ph/0701059v8 16 Jul 2007 A comparison of three turbulence models with an application to the West Pacific Warm Pool AC. BENNIS∗ ,M. GOMEZ MARMOL † ,R. LEWANDOWSKI ‡ ,T. CHACON REBOLLO §F. BROSSIER¶ Abstract In this work, we compare three turbulence models used to parameterize the oceanic boundary layer. These three models depend on the bulk Richardson number, which is coherent with the studied region, the West Pacific Warm Pool, because of the large mean shear associated with the equatorial undercurrent. One of these models, called R224, is new and the others are Pacanowski and Philander’s model (R213 model) and Gent’s model (R23 model). The numerical implementation is based on a non-conservative numerical scheme. The following (three criteria) are used to compare the models: the surface current intensity, the pycnocline’s form and the mixed layer depth. We initialize the code with realistic velocity and density profiles thanks the TOGA-TAO array (McPhaden, 1995, [21]). In case of static instability zone on the initial density profile, only the R224 model gives realistic results. Afterwards, we study a mixed layer induced by the wind stress. In this case, the R224 results and the Pacanowski and Philander’s results are similar. Furthermore, we simulate a long time case. We obtain a linear solution for all models that is in agreement with Bennis and al [1]. Summary 0.1 Keywords: vertical mixing, Richardson number, mixed layer. 1 Introduction The presence of an homogeneous layer near the surface of the ocean has been observed since a long time. This layer presents almost constant profiles of temperature and salinity. We distinguish the mixing layer from the mixed layer (Brainerd and Gregg, 1995, [3]). The mixing layer is actively mixed by surface fluxes. It is a depth zone where the turbulence is strong. The mixed layer is a maximum depth zone which has been mixed in the recent past by surface ∗IRMAR, Universit´e de Rennes 1, Campus de Beaulieu, 35042 Rennes Cedex, France †Departamento de Ecuaciones Diferenciales y An´alisis Numerico, Universidad de Sevilla. C/Tarfia, s/n.41080, Sevilla, Spain ‡IRMAR, Universit´e de Rennes 1, Campus de Beaulieu, 35042 Rennes Cedex, France §Departamento de Ecuaciones Diferenciales y An´alisis Numerico, Universidad de Sevilla. C/Tarfia, s/n.41080, Sevilla, Spain ¶IRMAR, Universit´e de Rennes 1, Campus de Beaulieu, 35042 Rennes Cedex, France 1 fluxes. This layer includes the mixing layer. In this work, we study the mixed layer. The bottom of the mixed layer corresponds either to the top of the thermocline, zone of large gradients of temperature, or to the top of the zone where haline stratification is observed (Vialard and Delecluse, 1998, [30]). Some attempts to describe this phenomenon can be found for example in Defant [5] or in Lewandowski [18]. The effect of the wind-stress acting on the sea-surface is considered to be the main forcing of this boundary layer. Observations in situ were completed by laboratory experiments (Deardoff, 1969, [4]) and more recently by numerical modelizations of the mixed layer. A historical of these different approaches can be found in K¨ase ([14]). K¨ase quoted the model of Kraus and Turner (1967, [15]) as the first applied mixed layer model which is a bulk one. There are also first order models as the Pacanowski and Philander model (called PP model, 1981, [26]) and the Large and Gent model (called KPP model, 1994, [16]). The second order models have been developed by Mellor and Yamada (called MY model, 1982, [22]) and Gaspar and al (1990, [7]). In this work, we focus our attention on the first order models. According to observations, the thickness of the mixed layer can vary between ten meters and a few hundred of meters, depending on the latitude (Boyer de Mont´egut and al, 2004, [23]). The mixed layer depth is difficult to determine. Traditionally, there are two types of criteria for this determination: the density difference criterion which will be used in this work and the density gradient criterion. The first one consists in estimating the mixed layer depth whose the density is about the surface density increased by 0.01 kg.m−3. This criterion is often used in the tropical area as in Peters and al, 1988 [27]. So, the mixed layer is the one where the density is inferior to the surface density to 0.01 kg.m−3. The second criterion is based on specifying density gradient. The mixed layer depth is where the density gradient is equal to the specified value. Mixing processes are intense in the homogeneous boundary layer, much weaker near the pycnocline or in the deep ocean. The effects of local fluxes of heat and momentum across the sea-surface induce turbulence and then mixing processes in the surface layer. The dynamical response of the ocean, and especially small scale turbulence, has a significant effect on mixing processes. This is true particularly in tropical areas. The rate of kinetic energy dissipation and the vertical turbulent fluxes of heat and mass cannot be measured directly, but they can be deduced from measurements of vertical temperature gradients and horizontal velocity [29]. Such experiments shown that turbulent dissipation is higher near the equator than in low latitudes [10], [11], [6]. Therefore in order to modelize the mixed boundary layer and to represent correctly the mixing processes, it is necessary to define a parametrization of turbulent diffusion. The mixed layer being strongly dominated by vertical fluxes, attention is drawn on vertical mixing which requires a closure model in order to represent the Reynolds stress. Recently, Goosse and al [9] studied the sensitivity of a global model to different parametrizations of vertical mixing. They insist on the crucial role of mixing in the upper oceanic layers: it has a direct impact on the sea surface temperature and then on ice evolution, but affects also the vertical profile of velocity by redistributing the momentum [2]. Classically, mixing parametrizations consist in the definition of turbulent eddy viscosity and 2 diffusivity coefficients. These coefficients can be chosen either as constants or as fixed profiles [13]. This is a simple but crude parametrization since variations of mixing with time and location are forbidden. A more suitable method is to define these coefficients as functions of processes governing the mixing. In tropical oceans vertical stratification and velocity shear are natural parameters since, following Philander [29], one of the reasons for higher turbulent dissipation near the equator than in low latitudes, is the large vertical shear observed in tropical currents. Pacanowski and Philander [26] proposed a formulation for eddy viscosity and diffusivity coefficients depending on the Richardson number which represents the ratio between buoyancy effects and vertical shear. The Richardson number dependent formulations allow strong mixing in high shear regions with low stratification, and low mixing elsewhere. It has to be noted that stratification tends to reduce turbulence and therefore mixing processes. In this paper, we study first order models. We focus on the behavior of Richardson number depending models: the basic one proposed by Pacanowski-Philander (PP model) and two variants used in Gent [8], Blanke and Delecluse [2],Goosse and al [9]. In this study, the PP model is called R213 model and the Gent model is called R23 model. The last model is named R224. The studied region is the West Pacific Warm Pool. At this location, the modelization based on the Richardson number is realistic because of the large mean shear associated with the equatorial undercurrent. Peters and al [27] have shown that the PP model underestimates the turbulent mixing at low Richardson number, while overestimating the turbulent mixing at high Richardson number in comparison with turbulent measurements. They found that the simulated thermocline is much too diffused compared with observations. Furthermore, this scheme overestimates the surface current intensity. Also, the simulated equatorial undercurrent is too shallow compared to observations ([12],[24],[19]). Halpern and al [12] compared the PP scheme and the MY scheme. At the equator, the current and temperature simulated by the PP scheme are more realistic than those simulated by the Mellor and Yamada scheme. The Gent model, one of three studied ones, gives realistic results in the West Pacific Warm Pool. This model simulates a sharp thermocline which is in agreement with the observations [8]. Furthermore, it gives good results for the annual average SST (Sea Surface Temperature) at the equator. In this paper, we investigate in case of each model, the simulated mixed layer depth, the form of the simulated pycnocline (sharp or diffuse) and the intensity of the simulated surface current. The numerical implementation is based on a non-conservative numerical scheme. We initialize the code with realistic velocity and density profiles. These profiles come from to the TOGATAO array (McPhaden, 1995, [21]). The geographic location is 0◦N, 165◦Ewhich is found in the West Pacific Warm Pool. In the first part, we study a mixed layer induced by the wind stress. Then, we initialize the code with a density profile showing a static instability zone. At last, we simulate a long time case. 3 2 Modelization of the mixed layer 2.1 Setting of the problem Variables used in this paper to describe the mixed layer are statistical means of the horizontal velocity and of the density denoted by (u, v, ρ).In the ocean, the density is a function of temperature and the salinity through a state equation. So, we consider the density as an idealized thermodynamic variable which is intended to represent temperature and salinity variations. The formation of the mixing layer is a response to sea-air interactions: wind-stress, solar heating, precipitations or evaporation acting on the sea surface. The variability of temperature is considered as essential to understand the response of the ocean. For example, the existence of a sharp pycnocline is a well known feature of the tropical areas. The thermal inertia of the water column is linked with the depth of the pycnocline, which influences the sea-surface temperature. The role of the haline stratification, sometimes considered as less important, has been recently evidenced by Vialard and Delecluse [30]: a numerical modelization of the tropical Pacific produces a ”barrier layer” depending on surface forcings and large scale circulation. The term ”barrier layer” refers to the water column located between the bottom of the mixed layer and the top of the pycnocline. It is present when the isohaline layer is shallower than the isothermal layer and then the depth of the mixed layer is controlled by the salinity. The model studied hereafter is not expected to describe all the phenomena occuring in the mixed layer. Its purpose is only a better understanding of a classical closure model. Therefore we shall use simplified equations governing the variables u(zonal velocity),v(meridian velocity) and ρ (density). The mixed layer being strongly dominated by vertical fluxes, we shall suppose that u,vand ρ are horizontally homogeneous and we so obtain a one-dimensional modelization. The Coriolis force will be neglected, which is valid in tropical oceans. Therefore, the equations governing the mixing layer are                ∂u ∂t =−∂ ∂z hu′w′i, ∂v ∂t =−∂ ∂z hv′w′i, ∂ρ ∂t =−∂ ∂z hρ′w′i, (1) where u′, v′, w′, ρ′represent the fluctuations of the horizontal velocity, vertical velocity and the density. The notation h i signifies that the quantity is statistically averaged. Equations (1) are the classical equations corresponding to a modelization of a water column. Such equations can be found in [17] as the equations of the boundary layer. Equations (1) are not closed and then the vertical fluxes appearing in the right-hand side have to be modelized. We study in this paper the behavior of a very classical closure modelization that uses the concept 4 of eddy coefficients in order to represent turbulent fluxes. So we set −u′w′=ν1 ∂u ∂z , −v′w′=ν1 ∂v ∂z , −ρ′w′=ν2 ∂ρ ∂z . Coefficients ν1and ν2are called vertical eddy viscosity and diffusivity coefficients and will be expressed as functions of the Richardson number Rdefined as R=−g ρ0 ∂ρ ∂z ∂u ∂z 2+∂v ∂z 2, where gis the gravitational acceleration and ρ0a reference density (for example ρ0= 1025 kg.m−3). The set of equations, initial and boundary conditions governing the mixing layer can now be written                                            ∂u ∂t −∂ ∂z ν1 ∂u ∂z = 0, ∂v ∂t −∂ ∂z ν1 ∂v ∂z = 0, ∂ρ ∂t −∂ ∂z ν2 ∂ρ ∂z = 0,for t⩾0 and −h⩽z⩽0, u=ub, v =vb, ρ =ρbat the depth z=−h, ν1 ∂u ∂z =ρa ρ0 Vx, ν1 ∂v ∂z =ρa ρ0 Vy, ν2 ∂ρ ∂z =Qat the surface z= 0, u=u0, v =v0, ρ =ρ0at initial time t= 0. (2) In system (2), the constant hdenotes the thickness of the studied layer that must contain the mixing layer. Therefore the circulation for z < −h, under the boundary layer, is supposed to be known, either by observations or by a deep circulation numerical model. This justifies the choice of Dirichlet boundary conditions at z=−h,ub,vband ρbbeing the values of horizontal velocity and density in the layer located below the mixed layer. The air-sea interactions are represented by the fluxes at the sea-surface: Vxand Vyare respectively the forcing exerced by the zonal windstress and the meridional wind-stress and Qrepresents the thermodynamical fluxes, heating or cooling, precipitations or evaporation. We have Vx=CD|ua|2and Vy=CD|va|2, where Ua= (ua, va) is the air velocity and CD(= 1,2.10−3) a friction coefficient. We study hereafter three different formulations for the eddy coefficients ν1=f1(R) and ν2= f2(R).Functions f1and f2can be defined as f1(R) = α1+β1 (1 + 5R)2, f2(R) = α2+f1(R) 1 + 5R=α2+α1 1 + 5R+β1 (1 + 5R)3.(3) 5 Formulation (3) corresponds to the modelization of the vertical mixing proposed by Pacanowski and Philander [26]. They proposed for coefficients α1, β1and α2the following values: α1= 1.10−4, β1= 1.10−2, α2= 1.10−5(units: m2s−1). This formulation has been used in the OPA code developed by Paris 6 University [2],[20] with coefficients α1= 1.10−6, β1= 1.10−2, α2= 1.10−7units: m2s−1.The selection criterion for the coefficients appearing in these formulas was the best agreement of numerical results with observations carried out in different tropical areas. A variant of formulation (3),proposed by Gent [8], is f1(R) = α1+β1 (1 + 10R)2, f2(R) = α2+β2 (1 + 10R)3(4) with α1= 1.10−4, β1= 1.10−1, α2= 1.10−5, β2= 1.10−1(units: m2s−1).A formulation similar to (4) when replacing 10Rby 5Rand varying the values of the coefficients α1, α2between the surface and the depth 50m is used in [9]. In this paper, we will also study the properties of another formulation close to formula (3): f1(R) = α1+β1 (1 + 5R)2, f2(R) = α2+f1(R) (1 + 5R)2=α2+α1 (1 + 5R)2+β1 (1 + 5R)4,(5) with α1= 1.10−4, β1= 1.10−2, α2= 1.10−5, β2= 1.10−3units: m2s−1. Eddy viscosity ν1defined by (5) is the same as the coefficient given by Pacanowski and Philander. The definition of the eddy diffusivity coefficient ν2differs by the exponent of the term (1 + 5R). In formulas (3) to (5), the eddy coefficients ν1and ν2are defined as functions of the Richardson number R through the terms (1 + γR)nappearing at the denominator. Hereafter, these three formulations will be denoted respectively by R213, R23 and R224 where R signifies Richardson number and the integer values are the exponents of (1 + γR) in the definitions of ν1and ν2. Eddy coefficients defined by relations (3) to (5) all present a singularity for a negative value of the Richardson number R=−0.2 or R=−0.1.We have plotted in Figure 1a the curves ν1=f1(R). In formulations (3) or (4) the coefficient of eddy diffusivity ν2becomes negative for values of Rlower than −0.2 or −0.1, and therefore the models are no more valid. The curves ν2=f2(R) obtained with formulations (3) (4) (5) and are plotted in Figure 1b. 6 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 Richardson Number Viscosity (m2/s) R=−0.2 R=−0.1 R224 R213 R23 Figure 1a: Viscosity (ν1=f1(R)) for all models −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −0.01 −0.008 −0.006 −0.004 −0.002 0 0.002 0.004 0.006 0.008 0.01 Richardson Number Diffusivity (m2/s) R=−0.2 R=−0.1 R224 R213 R23 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 Richardson Number Diffusivity (m2/s) R=−0.2 R=−0.1 R224 R213 R23 Figure 1b: Diffusivity (ν2=f2(R)) for all models Problem (2) coupled with one of the definitions (3) to (5) for eddy coefficients retains the vertical shear and buoyancy effect which are two important processes for the generation of the mixed boundary layer, especially in tropical areas. 7 3 Numerical Experiments 3.1 Finite Difference Scheme We want to resolve numerically the system (2). We replace the continuous variables (u, v, ρ, ν1, ν2) by discrete variables (un i, vn i, ρn i,(ν1)n i,(ν2)n i) which are the approximate solutions at time n∆t (with n= 1,2, ..., N) and at points (i−NI)∆z(with i= 1,2, ..., NI). We discretize the 1D-domain in z-levels where z is the vertical coordinate. We use a second-order central difference scheme for the second space derivative and a first-order backward difference scheme for the first space derivative. These previous schemes can be written as: •Second-order central difference scheme: ∂2u ∂z2n+1 i =un+1 i+1 −2un+1 i+un+1 i−1 ∆z2 •First-order backward difference scheme: ∂u ∂z n+1 i =un+1 i−un+1 i−1 ∆z The grid spacing, ∆z, is equal to 5 min cases 2 and 3 and equal to 1 min case 1. The time step, ∆t, is equal to 60 s. In time, we use implicit velocities and implicit density. The viscosity (ν1) and diffusivity (ν2) are explicit. The basin depth is 100 m. The boundary conditions are treated with a first-order backward difference scheme. We choose Neumann boundary conditions at the surface and Dirichlet boundary conditions in z=−h. The numerical scheme is the following:                                un+1 i−un i ∆t−ν1)n i−ν1)n i−1 ∆z· un+1 i−un+1 i−1 ∆z!−ν1)n i· un+1 i+1 −2un+1 i+un+1 i−1 ∆z2!= 0, vn+1 i−vn i ∆t−ν1)n i−ν1)n i−1 ∆z· vn+1 i−vn+1 i−1 ∆z!−ν1)n i· vn+1 i+1 −2vn+1 i+vn+1 i−1 ∆z2!= 0, ρn+1 i−ρn i ∆t−ν2)n i−ν2)n i−1 ∆z· ρn+1 i−ρn+1 i−1 ∆z!−ν2)n i· ρn+1 i+1 −2ρn+1 i+ρn+1 i−1 ∆z2!= 0. (6) and the Neumann boundary conditions at the surface are computed as: ν1)n i· un+1 i−un+1 i−1 ∆z!=ρa ρ0 Vx, ν1)n i· vn+1 i−vn+1 i−1 ∆z!=ρa ρ0 Vy, ν2)n i· ρn+1 i−ρn+1 i−1 ∆z!=Q. Moreover, the residual values are computed as: rn= NI X i=1 un+1 i−un i 2+ NI X i=1 vn+1 i−vn i 2+ NI X i=1 ρn+1 i−ρn i 2!1/2 .(7) 8 Our numerical scheme is non-conservative but it produces the same results as those given by a conservative one. Furthermore, the residual values decrease monotonically to zero as the flow reaches the steady states, better than the conservative scheme where the behaviour of the residual is more erratic. 3.2 Numerical Results In this section, we aim to study an Equatorial Pacific region called the West-Pacific Warm Pool. So, we initialize the code with data from the TAO array (McPhaden [21]). It is located at the equator between 120◦Eand 180◦E. The sea temperature is high and almost constant along the year (28 −30◦C). The precipitation are intense and hence the salinity is low. In the first part, we study a mixed layer induced by the wind stress. Then, we consider a case where the code is initialized by a density profile showing a static instability zone. Finally, we simulate a long time case. All our numerical results are grid-size (∆z) and time step (∆t) independent, in the sense that they remain practically unchanged as ∆zand ∆tdecrease. 3.2.1 The initial data We use data available from the Tropical Atmosphere Ocean (TAO) array (McPhaden [21]). The TAO project aims to study the exchange between the tropical oceans and the atmosphere. The TAO data have being very used in numerical simulations. The velocity data comes from the ACDP (Acoustic Doppler Current Profiler) measurements. To obtain the appropriate profiles, we interpolate the data by a one-order linear interpolation. We initialize our code by these profiles at 0◦N, 165◦Eand we obtain the results below. The buoyancy flux is equal to −1.10−6kg.m−2.s−1 (≃ −11 W/m2) in all cases. This heat flux is in agreement with Gent [8] because the heat flux must be in the range [0 W/m2; 20 W/m2] between 140◦E−180◦Eand 10◦N−10◦S. 3.2.2 Case 1: A mixed layer induced by the wind stress We use, in this case, as initial data, the TAO’s data for the time period between the 15th June 1991 and the 15th July 1991. The initial profiles are displayed on figure 2. The initial zonal velocity profile presents a westward current at the surface and, below it, an eastward undercurrent whose maximum is located about 55 m. Deepest, we observe a westward undercurrent. The initial density profile does not display a mixed layer. 9 and northward currents. The main southward current is located around to −55 mwhile the main northward current are located around to −90 mand at the surface. Furthermore, on initial density profile (see Figure 13), we notice a thirty five meters deep mixed layer according to Peters and al’s [28] criterion. −0.2 −0.15 −0.1 −0.05 0 0.05 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Initial zonal velocity (m/s) Depth (m) −0.12 −0.11 −0.1 −0.09 −0.08 −0.07 −0.06 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Initial meridian velocity (m/s) Depth (m) 1021.5 1022 1022.5 1023 1023.5 1024 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Initial density (kg/m3) Depth (m) Figure 13: The initial profile of zonal velocity, meridional velocity, density (left to right). The numerical results are displayed on figure 14. The three turbulence models give a linear profile for the simulated time. This fact corroborates the existence of a linear equilibrium solution obtained by Bennis and al [1]. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Final zonal velocity (m/s) Depth (m) R224 R213 R23 −0.105 −0.1 −0.095 −0.09 −0.085 −0.08 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Final meridian velocity (m/s) Depth (m) R224 R213 R23 1023.495 1023.5 1023.505 1023.51 1023.515 1023.52 1023.525 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Final density (kg/m3) Depth (m) R224 R213 R23 Figure 14: Comparison of zonal velocity profiles (left position), meridional velocity profiles (medium position) and density profiles (right position). The residual values are displayed on figure 15. Notice that the numerical convergence to steady states is good. 16 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 1 2 3 4 5 6x 10−4 Time (h) Residual R224 R213 R23 Figure 15: Comparison of residual values. As the equilibrium solutions (zonal velocity ue, meridian velocity veand density ρe) are unique in case of the R224 model (see Bennis and al, 2007 [1]), we can compare the theoretical solutions with the numerical solutions. We remember hereafter the theoretical equilibrium solutions which satisfy the following system: ∂ ∂z f1(Re)∂ue ∂z = 0,∂ ∂z f1(Re)∂ve ∂z = 0,∂ ∂z f2(Re)∂ρe ∂z = 0. The theoretical steady-state profiles for velocities and density are obtained by integrating the previous system with respect to z, taking account the boundary conditions at z=−h: ue(z) = ub+Vxρa ρ0f1(Re)(z+h), ve(z) = vb+Vyρa ρ0f1(Re)(z+h), ρe(z) = ρb+Q f2(Re)(z+h). The equilibrium richardson number (Re) can be interpreted as the intersection of the curves k(R) = (f1(R))2 f2(R)and h(R) = CR with C=−ρ2 a(V2 x+V2 y) gQρ0 . In case of R224 model, the graph of function kand hfor Q < 0 is plotted on figure 16. −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 R224 Richardson Number Re =0.063935 h(R) k(R) Figure 16: Equilibrium richardson number - R224 model. 17 So, we compute the theoretical solutions with the equilibrium richardson number (Re = 0.063935) given by the equation h(R) = k(R). The theoretical and numerical equilibrium solutions are displayed on figure 17. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Final zonal velocity (m/s) Depth (m) R224 Theoretical Solution Numerical Solution −0.105 −0.1 −0.095 −0.09 −0.085 −0.08 −0.075 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Final meridian velocity (m/s) Depth (m) R224 Theoretical Solution Numerical Solution 1023.495 1023.5 1023.505 1023.51 1023.515 1023.52 1023.525 −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 Final density (kg/m3) Depth (m) R224 Theoretical Solution Numerical Solution Figure 17: Theoretical and numerical equilibrium solutions - R224 model. We obtain similar profiles and the small differences are due to numerical scheme which we use. 4 Conclusion 4.1 Code validation We observe, in case 1, a mixed layer induced by the wind stress that is a physical phenomenon. Moreover, we notice, in all cases, an increase in the zonal surface current when we apply an eastward wind (u > 0) at the surface. In the same order, an efficiency northward wind (v > 0) at the surface causes an increase in the meridional surface current. These results are in agreement with the physical reality. Furthermore, the residual values correctly converge to steady states, corresponding to rather high levels of turbulent diffusion. Theses previous observations validate our code that gives a good representation of the physical reality. 4.2 Comparison of three turbulence models Our comparison is based on three criteria: the mixed layer depth, the surface current intensity and the pycnocline’s form. The mixed layer depth is obtained with a density difference criterion. This difference is equal to 0.01 kg.m−3. We use the surface current values to determine the surface current intensity. The density gradient is used to determine the pycnocline’s form. We summarize the results hereafter. - In the particular case 1, the R213, R23 and R224 models give a same mixed layer depth. In terms of surface current, those simulated by R224 and R213 models are of the same kind. However the R224 surface current is slightly stronger than the R213 surface current. The 18 R23 model underestimates this current. Furthermore, the pycnocline’s form are similar for R213, R23 and R224 models. - In the particular case 2, the R213 and R23 models have negative diffusivities at the initial time. Therefore, these models are not physically valid in this case. This problem comes from static instability zone in the initial density profile. Hence, we can not use these models in this case. Then, only R224 is valid. The R224 model produces a homogeneous mixed layer. So, we do not compare the pycnocline’s form since we have only one mixed layer. - Furthermore, the long time profiles are linear for all models that is in agreement with Bennis and al [1]. In case of R224 model, the theoretical and the numerical profiles are similar. Globally, the R224 model has the same behavior as the Pacanowski and Philander model (R213 model) and we can use it in more situation. For example, R224 model can be used in almost all of convective cases. Moreover, the R224 numerical equilibrium solution is in agreement with the R224 theoretical equilibrium solution. Acknowlegments The research of T. Chacon Rebollo and Macarena Gomez Marmol has been partially funded by Spanish Government Grant MTM2006-01275. Furthermore, the authors express their gratitude to Pascale Delecluse, Assistant Director of Research at Meteo-France, for her relevant advice. References [1] A. C. Bennis, T. C. Rebollo, M. G. Marmol, and R. Lewandowski,Stability of some turbulent vertical models for the ocean mixing boundary layer, Applied Mathematical Letters, To Appear (2007). [2] B. Blanke and P. Delecluse,Variability of the tropical atlantic ocean simulated by a general circulation model with two different mixed-layer physics, J. Phys. Oceanography, 23 (1993), pp. 1363–1388. [3] K. Brainerd and M. Gregg,Surface mixed layer and mixing layer depths, Deep Sea research, 42 (1995), pp. 1521–1543. [4] J. W. Deardorff, G. Willis, and D. Lilly,Laboratory investigation of nonsteady penetrative convection, J. Fluid. Mech., 35 (1969), pp. 7–31. [5] A. Defant,Schichtung und zirculation des atlantischen ozeans, Wiss. Ergebn. Deutsch. Atlant. Exp. Meteor, 6 (1936), pp. 289–411. [6] A. Gargett and T. R. Osborn,Small scale shear measurements during the fine and microstructure experiment, J. Geophys. Res., 86 (1981), pp. 1929–1944. 19 [7] P. Gaspar, Y. Gregoris, and L. J. M.,A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: test at station papa and long-term upper ocean study site, J. Geophys. Research, 16 (1990), pp. 179–193. [8] P. R. Gent,The heat budget of the toga-coare domain in an ocean model, J. Geophys. Res., 96 (1991), pp. 3323–3330. [9] H. Goosse, E. Deleersnijder, T. Fichefet, and M. H. England,Sensitivity of a global coupled ocean-sea ice model to the parametrization of vertical mixing, J. Geophys. Res., 104 (1999), pp. 13681–13695. [10] M. C. Gregg,Temperature and salinity microstructure in the pacific equatorial undercurrent, J. Geophys. Res., 81 (1976), pp. 1180–1196. [11] M. C. Gregg and T. B. Sandford,Signature of mixing from the bermuda slope, the sargasso sea and the gulf stream, J. Phys. Oceanogr., 10 (1980), pp. 105–127. [12] D. Halpern, Y. Chao, C. Ma, and C. Mechoso,Comparison of tropical pacific temperature and current simulations with two vertical mixing schemes embedded in ocean general circulation model and reference to observations, J. Geophys. Res, 100 (1995), pp. 2515–2523. [13] A. C. Hirst and T. J. M. Dougall,Deep-water properties and surface buoyancy flux as simulated by a z-coordinate model including eddy-induced advection, J. Phys. Oceanogr., 26 (1996), pp. 1320–1343. [14] R. H. K¨ ase,Modeling of the oceanic mixed-layer and effects of deep convection, 1998. [15] E. Kraus and J. Turner,A one dimensional model of the seasonal thermocline, ii, the general theory and its consequences, Tellus, 19 (1967), pp. 99–105. [16] W. G. Large, C. McWilliams, and S. C. Doney,Oceanic vertical mixing : a review and a model with a nonlocal boundary layer parametrization, Rev. Geophys., 32 (1994), pp. 363–403. [17] M. lesieur,Turbulence in Fluids, Kluwer, 1997. [18] R. Lewandowski,Analyse math´ematique et oc´eanographie, Masson, 1997. [19] X. Li, Y. Chao, J. McWilliams, and L. Fu,A comparison of two vertical-mixing schemes in a pacific ocean general circulation model, Journal of Climate, 14 (2001), pp. 1377–1398. [20] G. Madec, P. Delecluse, M. Imbard, and C. Levy,O.p.a. version 8.0. ocean general circulation model, reference manual, 1997. Technical report. [21] M. McPhaden,The tropical atmosphere ocean (tao) array is completed, Bull. Am. Meteorol. Soc, 76 (1995), pp. 739–741. 20 [22] G. Mellor and T. Yamada,Development of a turbulence closure model for geophysical fluid problems, Reviews of Geophysics and Space Physics, 20 (1982), pp. 851–875. [23] C. B. Montegut, G. Madec, A. S. Fischer, A. Lazar, and D. Iudicone,Mixed layer depth over the global ocean : An examination of profile data and a profile-based climatology, Journal Geophysical Research, 109 (2004), p. C12003. [24] P. P. Niiler and Coauthors,Comparison of toga tropical pacific ocean model simulations with the woce/toga surface velocity programme drifter data set., World Climate Research Programme Rep. WCRP-1995, (1995), p. 156pp. [25] T. Osborn and C. Cox,Oceanic finestructure, Geophys. Fluid Dyn., 3 (1972), pp. 321– 345. [26] R. C. Pacanowski and S. G. H. Philander,Parametrization of vertical mixing in numericals models of the tropical oceans, J. Phys. Oceanogr., 11 (1981), pp. 1443–1451. [27] H. Peters, M. C. Gregg, and J. M. Toole,On the parametrization of equatorial turbulence, Journal of Geophysical Research, 93 (1988), pp. 1199–1211. [28] ,Meridional variability of turbulence through the equatorial undercurrent, Journal of Geophysical Research, 94 (1989), pp. 18,003–18,009. [29] S. G. H. Philander,La Ni˜na, and the Southern Oscillation, Academic Press, 1990. [30] J. Vialard and P. Delecluse,An ogcm study for the toga decade. part i: Role of salinity in the physics of the western pacific fresh pool, J. Phys. Oceanogr., 28 (1998), pp. 1071–1088. 21