Abstract

The formation of a diurnal thermocline in the ocean mixed layer under a stabilizing buoyancy flux was simulated successfully by large-eddy simulation, reproducing various features consistent with observation. The analysis of the simulation result revealed that the formation of a diurnal thermocline passes through two different phases: the formation of a thermocline (formation stage) and increasing thickness of the thermocline thereafter (growth stage). Turbulent kinetic energy (TKE) flux dominates TKE production within the mixed layer, but turbulence maintained by shear production at the thermocline causes stratification below the mixed layer. In addition, once the thermocline is formed, both the gradient and flux Richardson numbers maintain constant values at the thermocline. It was also found that a diurnal thermocline cannot be formed in the absence of both wave breaking and Langmuir circulation. Furthermore, the effects of stratification on turbulence were investigated based on the time series of various physical variables of turbulence at the diurnal thermocline and within the mixed layer, and the mechanism for diurnal thermocline formation is discussed.

1. Introduction

Strong turbulence usually exists near the surface in the ocean mixed layer as a result of wave breaking (WB; e.g., Agrawal et al. 1992; Drennan et al. 1996), leading to the response to a surface stabilizing buoyancy flux that is fundamentally different from the atmospheric boundary layer. A diurnal thermocline (or “thermocline” hereafter) is formed at a certain depth during the day in the ocean mixed layer while a temperature gradient remains small near the surface. A strong temperature gradient, however, appears near the surface during the night in the atmospheric boundary layer.

Understanding of the dynamical process of diurnal thermocline formation is essential for predicting the sea surface temperature and the vertical transport of heat, momentum, and dissolved gases in the upper ocean. Only a few studies, however, have been made so far to understand the diurnal thermocline, based on field observation (Stommel et al. 1969; Delnore 1972; Kudryavtsev and Soloviev 1990; Brainerd and Gregg 1993; Caldwell et al. 1997), mixed layer models (Woods and Barkmann 1986; Kraus 1988; Noh and Fernando 1991; Noh 1996), or laboratory experiments (Kantha and Long 1980; Hopfinger and Linden 1982; Noh and Long 1990), and the mechanism for its formation is not yet clearly understood.

Noh (1996) suggested that turbulent kinetic energy (TKE) production is dominated by the TKE flux divergence F in the upper mixed layer as a result of WB, and that a diurnal thermocline is formed by the positive feedback between F and buoyancy decay Pb. He also suggested that once the thermocline is formed, local balance is observed among Pb, shear production Ps, and dissipation ɛ (i.e., PsPbε = 0), and the flux Richardson number Rf (= Pb/Ps) remains constant at the thermocline. Turbulence maintained by shear production at the thermocline causes stratification below the mixed layer by allowing heat transport across the thermocline, which is also found in observation data (Brainerd and Gregg 1993; Caldwell et al. 1997). These suggestions were based on the mixed layer model results, however.

Meanwhile, recent progress in the reproduction of realistic turbulent flows in the ocean mixed layer by using large-eddy simulation (LES: Skyllingstad and Denbo 1995; McWilliams et al. 1997; Skyllingstad et al. 1999; Noh et al. 2004, 2006; Min and Noh 2004; Li et al. 2005; Sullivan et al. 2007; Polton and Belcher 2007) provides us with the possibility of investigating the dynamical process of the ocean mixed layer at a fundamental level. Although there are a few examples of LES for the ocean mixed layer under a stabilizing buoyancy flux (Min and Noh 2004; Li et al. 2005), the formation of a diurnal thermocline has not been considered yet.

Therefore, in this paper an attempt has been made to reproduce the formation of a diurnal thermocline in the ocean mixed layer under a stabilizing buoyancy flux by using LES and to clarify its dynamical process by analysis of LES data. In particular, investigation was focused on how turbulence at the thermocline is modified during the formation of a diurnal thermocline. It was also investigated how the result is affected by Langmuir circulation (LC), WB, radiation penetration, and the diurnal variation of the surface buoyancy flux.

2. Simulation

For the simulation, we used the LES model for the ocean mixed layer developed by Noh et al. (2004), in which both LC and WB are realized. The modified filtered equation is given by

 
formula

where π is the generalized pressure [= p/ρ0 + (|u + us|2 − |u|2)/2], τij is the subgrid scale Reynolds stress, and b is buoyancy [= −g(ρρ0)/ρ0]. Here LC is generated by the vortex force ɛijkusjωk with the Stokes drift velocity usi (Craik and Leibovich 1976), and WB is realized by a stochastic forcing Xi.

The random forcing Xi is designed to reproduce the turbulence with the length l0 and time scale ℑ0 corresponding to those of near-surface small-scale turbulence generated by WB (Noh et al. 2004). It is expressed as

 
formula

where G(0;1) is the Gaussian random function with a mean of 0 and variance of 1, δ is the delta function, and u* is the frictional velocity. The length and time scales of random forcing at the surface are given by l0 = 1.25 m and ℑ0 = 0.1 l0/Au*, based on observational evidence (Agrawal et al. 1992; Craig and Banner 1994; Drennan et al. 1996; Gemmrich and Farmer 1999). The value of an empirical constant A = 3 is determined so that the profile of the resultant dissipation rate is consistent with observation. Here the random forcing is designed to reproduce the realistic turbulence structure near the sea surface generated by WB rather than to simulate a realistic WB process.

For simplicity, we assumed that both the wind stress and wave fields are in the x direction and further assumed that the wave field is steady and monochromatic. The associated Stokes velocity is then given by us = Us exp(−4πz/λ), with Us = (2πa/λ)2(/2π)1/2, where a is the wave height, λ is the wavelength, and g is the gravitational acceleration. For wave height and wavelength, typical values were used such as a = 1.0 m and λ = 40 m, which makes Us = 0.196 m s−1. Simulations were also performed in the absence of WB (A = 0) and/or LC (Us = 0) to investigate their effects.

The LES model used in this study was developed based on the Parallelized LES Model (PALM), which has been applied to simulate various geophysical turbulence phenomena. Details of the LES code can be found in Raasch and Schröter (2001). The model domain was 300 m in the horizontal direction (x and y) and 80 m in the vertical direction (z). The number of grid points was 240 × 240 × 64, and the corresponding grid sizes were 1.25 m in all directions. Simulation was also performed with 240 × 240 × 128 grid points (Δz = 0.625 m) to examine the sensitivity to the vertical resolution. A free-slip boundary condition was applied at the bottom. The wind stress was given by u* = 0.01 m s−1, which corresponds to the wind velocity 5–10 m s−1 at 10 m above the sea level (Fairall et al. 2003). The Coriolis force was given by f = 0.5 × 10−4 s−1, corresponding to latitude 20°, so the Ekman length scale u*/f is much larger than the Monin–Obukhov length scale u*3/Q0.

Initially, integration started with a mixed layer of uniform density without surface buoyancy flux. Once a quasi-equilibrium state was reached after an 8-h integration, the constant surface buoyancy flux Q0 = 5 × 10−7 m2 s−3 was imposed. Turbulence reaches quasi-equilibrium over a period of O( f−1) in the LES of the ocean mixed layer (McWilliams et al. 1997). This moment was defined as the initial time (i.e., t = 0 h). If the surface freshwater flux and salinity variation are neglected, Q0 corresponds to the surface heat flux H0 by Q0 = (/ρcp)H0, where α = −ρ−1/dT and cp is the specific heat at constant pressure, and a pycnocline corresponds to a thermocline. The condition in the present simulation corresponds to Lat [≡ (u*/Us)1/2] = 0.23 and Ho [≡ −2Q0λ/(2πUsu*2)] = −0.33.

Because the primary concern in this paper is to understand the fluid dynamical process of diurnal thermocline formation, a constant buoyancy flux was used and solar radiation penetration was neglected. The effects of radiation penetration and diurnal variation of the surface heat flux were investigated, however, using the simulations considering these effects.

3. Results

a. Formation of a diurnal thermocline

Figure 1 shows how buoyancy distribution in the vertical cross section evolves in the initial stage of thermocline formation. At t = 0.5 h, a well-defined thermocline (or pycnocline), has not formed yet, and buoyancy is transported downward (similar to the case of a passive scalar). Buoyancy penetrates deeper at some locations, following downward jets of LC. The front of buoyancy discontinuity appears at t = 1 h at a certain depth, where fluctuation continues to decrease with increasing stratification.

Fig. 1.

Variation of buoyancy distribution at the vertical cross section with time (Δt = 0.5 h).

Fig. 1.

Variation of buoyancy distribution at the vertical cross section with time (Δt = 0.5 h).

Corresponding evolutions of the horizontal mean values of buoyancy B, vertical shear of the horizontal mean velocity S2 [= (∂U/∂z)2 + (∂V/∂z)2], and dissipation rate ɛ are shown in Fig. 2, but for a longer period (t < 8 h). Hereafter all variables appearing in the present paper are horizontally averaged ones, and dissipation and flux terms include both resolved and subgrid-scale parts.

Fig. 2.

Evolution of vertical profiles with time (Δt = 1 h, and the dotted lines represent t = 0 and 0.5 h): (a) B, (b) S2, and (c) ɛ.

Fig. 2.

Evolution of vertical profiles with time (Δt = 1 h, and the dotted lines represent t = 0 and 0.5 h): (a) B, (b) S2, and (c) ɛ.

An exponentially decreasing B appears at early times (t = 0.5 h), similar to the case of constant eddy diffusivity. After some time (t = 1 h), a weak thermocline appears near z = 8 m and stratification N2 (= ∂B/∂z) at this depth continues to increase until t = 2 h. This depth will be called the mixed layer depth h in the present paper. After t = 2 h, however, the thermocline increases its thickness with time, causing stratification below the mixed layer, whereas the increasing rate of N2 becomes much weaker (see also Fig. 4a).

Fig. 4.

Time series of physical variables (dotted line: z = 3.75 m; solid line: z = 7.5 m; dashed line: z = 12.5 m): (a) N2, (b) S2, (c) Rf, (d) Ri, (e) /2, (f) ( + )/2, (g) /( + ), (h) Kh, (i) Km, (j) lm, (k) lb, (l) ld, and (m) ls.

Fig. 4.

Time series of physical variables (dotted line: z = 3.75 m; solid line: z = 7.5 m; dashed line: z = 12.5 m): (a) N2, (b) S2, (c) Rf, (d) Ri, (e) /2, (f) ( + )/2, (g) /( + ), (h) Kh, (i) Km, (j) lm, (k) lb, (l) ld, and (m) ls.

From the balance of the TKE budget, Niiler and Kraus (1977) predicted that h is proportional to the Monin–Obukhov length scale L (≡u*3/Q0) or h = mL with the proportional constant m in the range 2.5–16, although they did not predict the growth of thermocline thickness. It gives h ∼ 5–32 m in the present simulation, which is consistent with Fig. 2a.

Formation of a diurnal thermocline suppresses the momentum flux as well as the buoyancy flux across it, and it induces velocity shear at the thermocline. The vertical shear S2 appears at t = 1 h and increases with time. After t = 2 h, however, the vertical range of S2, which is equivalent to the thickness of the thermocline, increases with time, whereas the increasing rate of S2 for a given depth becomes much weaker (see also Fig. 4b). It is also observed that increasingly larger values of the maximum S2 appear at deeper depth with time (see also Fig. 4b).

The dissipation rate ɛ decreases below the thermocline initially, as the TKE flux from the sea surface is shut off by the formation of a thermocline, and then increases later at the thermocline (10 m < z < 20 m) after t > 2 h along with increasing shear production.

The evolution of profiles of B, S, and ɛ (shown above) suggests that the formation of a diurnal thermocline passes through two different phases: the formation of a thermocline until t ∼ 2 h (“formation stage”) and increasing thickness of the thermocline thereafter (“growth stage”).

Most features shown in Fig. 2 are clearly evidenced from observational data (Brainerd and Gregg 1993; Caldwell et al. 1997). For example, Fig. 3 shows the formation of a diurnal thermocline at z ∼ 10 m and the decrease of ɛ below the mixed layer with ɛ/(dɛ/dt) ∼1–2 h during the formation stage (1228 LT), followed by the growth of thermocline thickness up to 10 m and the slow increase of ɛ at the thermocline during the growth stage (1423 LT, 1646 LT). Note that the formation of the thermocline is slower in Fig. 3 because of the diurnal variation of Q0 starting with a very small value (see section 3f). The profiles equivalent to Fig. 2 could also be obtained from the mixed layer model by Noh and Kim (1999, their Fig. 2).

Fig. 3.

Profiles of hourly averages of dissipation ɛ (shaded) and potential temperature θ on 17 Oct during Patches Experiment (PATCHEX) (Brainerd and Gregg 1993, their Fig. 12). The surface stabilizing buoyancy flux was applied during 0900 and 1623 LT. (1 MPa = 100 m in water)

Fig. 3.

Profiles of hourly averages of dissipation ɛ (shaded) and potential temperature θ on 17 Oct during Patches Experiment (PATCHEX) (Brainerd and Gregg 1993, their Fig. 12). The surface stabilizing buoyancy flux was applied during 0900 and 1623 LT. (1 MPa = 100 m in water)

b. Evolution of turbulence at the diurnal thermocline

To investigate how turbulence is modified during the formation of the diurnal thermocline, the time series of various physical variables were calculated at z = 3.75, 7.5, and 12.5 m (Fig. 4). Here z = 3.75 m represents the depth within the mixed layer. The thermocline starts to form near z = 7.5 m and reaches z = 12.5 m at a later time as its thickness increases. Time series analysis confirms the transition from the formation stage to the growth stage at t ∼ 2 h.

At z = 7.5 m, N2 increases rapidly initially, but the increasing rate becomes much smaller after t ∼ 2 h. The development of N2 is delayed at deeper depth (z = 12.5 m). It is also found that N2 approaches a larger value at deeper depth, in the same way as S2 (Fig. 2b), which is due to weaker TKE and vertical mixing there (see Figs. 4e,h). The development of S2 follows a similar pattern but lags behind N2 by 0.5–1 h. Note that S2 (z = 7.5 m) is still negligible at t = 0.5 h. Both N2 and S2 remain very small within the mixed layer (z = 3.75 m), as expected from Fig. 2, although they also increase slowly with time. The values of N2 and S2 are comparable to observational data both within the mixed layer and at the thermocline (Brainerd and Gregg 1993; Caldwell et al. 1997).

Both Rf (= Pb/Ps) and Ri [= (N/S)2] at the thermocline increase initially to very large values (Rf ≅ 2 and Ri ≅ 4 near t = 0.5 h) and then fall to constant values of Rf ≅ 0.4 and Ri ≅ 0.35. The initial peaks of Rf and Ri can be understood by the time lag in the growth of S2 compared to N2, as shown in Figs. 4a,b. Caldwell et al. (1997) observed Ri ∼ 0.2–0.5 at the diurnal thermocline. Both Ri and Rf are always larger than 1 within the mixed layer in which Ps is negligible (see also section 3e).

Both horizontal and vertical TKE, /2 and ( + )/2, at the thermocline decrease rapidly during the formation stage but increase slowly thereafter. The vertical TKE is more strongly affected by stratification than the horizontal TKE so that the ratio r [= /( + )] drops rapidly to about 0.15 during the initial stage, but it remains invariant thereafter. It is consistent with observation data r ∼ 0.15–0.25 (Schumann and Gerz 1995; Caughey et al. 1979). Within the mixed layer, both vertical and horizontal TKE show a weak decreasing tendency, and r remains at a much larger value of r ≅ 0.32.

Vertical eddy viscosity Km and diffusivity Kh are calculated by (2 + 2)1/2 = KmS and − = KhN. Their values at the thermocline become very large initially because momentum and buoyancy are effectively transferred downward by large-scale eddies such as LC in this stage, even though N2 and S2 are negligible. After the formation of a thermocline, Km and Kh fall by a factor of 10–100. Much stronger vertical mixing is always maintained within the mixed layer.

Finally, we examined the mixing length scale lm, defined by Kh = qlm, and the dissipation length scale ld, defined by ɛ = q3/ld, using q = ()1/2, in comparison with lb (= q/N) and ls (= q/S), which are often used for the scaling of length scales (Figs. 4j–m; Britter et al. 1983; Hunt et al. 1988). Although lm shows a similar temporal variation to lb for a given depth, lm/lb is different depending on depth. It suggests that lm cannot be estimated by q and N only but may also be influenced by other factors, such as z.

On the other hand, ld reveals a temporal variation that is apparently unrelated to lb. The initial decrease of ld at z > 7.5 m can be understood from the fact that large eddies from the surface, including LC, are broken down by the formation of a thermocline rather than scaled by lb or ls.

c. TKE budget

The TKE budget in the presence of LC and WB can be derived from (1) as

 
formula

under the horizontally homogeneous condition, where E (= /2) is TKE, F [= −∂Π/∂z = −∂( + /2 − 2υ)/∂z] is the divergence of TKE flux Π, Ps (= −U/∂zV/∂z) is shear production, Pb (= −) is decay by buoyancy, ɛ is dissipation, PL [= us∂()/∂z] is the contribution from the vortex force (see, e.g., Kantha and Clayson 2000), and W is the TKE influx at the surface from wave breaking. Here ui, π′, and b′ are fluctuating components of velocity, generalized pressure, and buoyancy, and sij = (∂ui/∂xj + ∂uj/∂xi)/2. Noh et al. (2004) showed that W from a random forcing Xi in (1) is evaluated as αu*3 with α ≅ 40. Because W dominates TKE production near the surface, the integration of (3) over an arbitrarily small depth δ from the surface leads to Π(δ) ≅ W from Π(0) = 0. Therefore, W can be regarded as the surface boundary condition of Π rather than a separate term in the TKE budget, if we only consider the TKE budget below z = δ.

Profiles of the terms in the TKE budget in the ocean mixed layer at t = 8 h, shown in Fig. 5, reveal that TKE production is dominated by F within the mixed layer (z < h), whereas the local balance PsPbε = 0 is observed at the thermocline, as suggested by Noh (1996). It is also found that PL is much smaller than F within the mixed layer and becomes negligible at the thermocline.

Fig. 5.

Profiles of the terms of TKE budget (F, thick solid; Ps, dashed; Pb, dotted; ɛ, thin solid; PL, dash–dot–dash).

Fig. 5.

Profiles of the terms of TKE budget (F, thick solid; Ps, dashed; Pb, dotted; ɛ, thin solid; PL, dash–dot–dash).

d. Effects of wave breaking and Langmuir circulation

Simulations were also performed in the absence of WB and/or LC to investigate their effects. Figure 6 shows the evolution of B profiles from experiments without both LC and WB (EXP O), with WB only (EXP W), and with LC only (EXP L).

Fig. 6.

Evolution of vertical profiles of B with time (Δt = 1 h, and the dotted line represents t = 0.5 h): (a) EXP O, (b) EXP W, and (c) EXP L.

Fig. 6.

Evolution of vertical profiles of B with time (Δt = 1 h, and the dotted line represents t = 0.5 h): (a) EXP O, (b) EXP W, and (c) EXP L.

Strong stratification appears near the surface without forming a diurnal thermocline in EXP O, similar to the atmospheric boundary layer (Fig. 6a). In this case, the maximum values of both S2 and N2 appear near the surface. It is important, however, to mention that the case corresponding to EXP O actually occurs in the ocean under the very weak wind stress in which WB and LC cannot occur. For example, Soloviev and Lukas (1997) observed that strong stratification near the sea surface appears during calm weather in the western equatorial Pacific warm pool. A diurnal thermocline is generated in EXP L and EXP W (Figs. 6b,c), although the depth of the thermocline from EXP W and EXP L is slightly shallower than in the case where both WB and LC are present (EXP LW; Fig. 2a). It should be mentioned, however, that, if Ho becomes very small (Ho < −0.5 at Lat = 0.23), LC breaks down and strong stratification appears near the surface in the case of LC only, similar to the case of EXP O (Min and Noh 2004).

The corresponding profiles of the terms of the TKE budget are shown in Fig. 7. In EXP O, shear production dominates TKE production, as expected, as in the atmospheric boundary layer, and PsPbε = 0 is observed from the surface. The case from EXP W is similar to that from EXP LW, although shear production is somewhat larger at the thermocline. On the other hand, EXP L suggests that the vertical transport of buoyancy and TKE by LC may play a similar role to the TKE flux. The existence of a small amount of F near the surface in EXP L reflects the strong TKE production near the surface by PL, which can be inferred from the substantial increase of vertical TKE near the surface in the presence of LC (McWilliams et al. 1997; Noh et al. 2004; Li et al. 2005) and its downward transport by LC.

Fig. 7.

Profiles of the terms of TKE budget: (a) EXP O, (b) EXP W, and (c) EXP L (F, thick solid; Ps, dashed; Pb, dotted; ɛ, thin solid; PL overlaps with F in EXP L).

Fig. 7.

Profiles of the terms of TKE budget: (a) EXP O, (b) EXP W, and (c) EXP L (F, thick solid; Ps, dashed; Pb, dotted; ɛ, thin solid; PL overlaps with F in EXP L).

The comparison of the time series of Rf, Ri, and r at the thermocline (z = 12.5 m) shows that the large values of these parameters at the initial stage appear only in EXP L and EXP LW (Fig. 8). It explains that these large values are due to LC, which causes the strong vertical component of TKE and strong vertical mixing without shear. On the other hand, the equilibrium values remain invariant, suggesting that they are determined by the interaction between turbulence and stratification and insensitive to the details of WB and LC.

Fig. 8.

Time series of physical variables at z = 12.5 m (EXP LW, thick solid; EXP L, dashed; EXP W, thin solid; EXP O, dotted): (a) Rf, (b) Ri, and (c) /( + ).

Fig. 8.

Time series of physical variables at z = 12.5 m (EXP LW, thick solid; EXP L, dashed; EXP W, thin solid; EXP O, dotted): (a) Rf, (b) Ri, and (c) /( + ).

e. Mechanism for the formation of a diurnal thermocline

Important questions arising from the diurnal thermocline formation shown in previous sections are why a thermocline appears at a certain depth (formation stage) and why the thermocline thickness grows with time thereafter (growth stage).

The fact that F disappears at z = h (Fig. 5) suggests the role of F in the formation of a diurnal thermocline, whereas the growth of the thermocline thickness occurs under the condition PbPs − ɛ = 0. Based on these features, each stage can be explained in terms of a different feedback mechanism, as summarized in Table 1 (Noh 1996).

Table 1.

Feedback mechanisms between Pb vs F and Pb vs Ps.

Feedback mechanisms between Pb vs F and Pb vs Ps.
Feedback mechanisms between Pb vs F and Pb vs Ps.

According to the feedback mechanism presented in Table 1, under the interaction between Pb and F, TKE and Kh decrease and N increases continuously with time at a certain depth, leading to the formation of a diurnal thermocline across which fluxes of both buoyancy and TKE are prohibited. On the other hand, under the interaction between Pb and Ps, the values of TKE, Kh and N are maintained at certain levels, and it allows the downward transport of buoyancy below the mixed layer.

These two different feedback mechanisms are also reflected in the evolution of ɛ shown in Fig. 2c. During the formation stage (t < 2 h), the positive feedback between F and Pb causes the decrease of ɛ with time, but during the growth stage (t > 2 h) it increases while approaching PsPb − ɛ = 0 along with the increase of Ps. Noh (1996) also showed using the mixed layer model that the thermocline cannot increase its thickness in the absence of Ps.

The relation Rf = 1 − ɛ/Ps can be obtained from PsPb − ɛ = 0 at the thermocline. When Ps is the only TKE source, ɛ = msPs with ms ∼ 0.5 was suggested by Niiler and Kraus (1977) and Davis et al. (1981), which is also supported from Fig. 5. It implies Rf ∼ 0.5 at the thermocline. On the other hand, within the mixed layer, the relation Rf = 1 − ɛ/Ps + F/Ps is obtained, and a much larger value of Rf is expected because FPs there. Figure 4c shows that Rf (z = 3.25 m) ≅ 1.2, Rf (z = 7.5 m) ≅ 0.5, and Rf (z = 7.5 m) ≅ 0.4. Because Pr (≡ Km/Kh) remains invariant with time (Figs. 4h,i), Ri [= Rf(Pr)] also remains constant (Fig. 4d).

f. Effects of radiation penetration and diurnal variation of heat flux

In the present work, the effects of radiation penetration have been neglected because it is not directly involved in the mechanism for the formation of diurnal thermocline by the interaction between turbulence and stratification. Nonetheless, the diurnal thermocline is usually formed at a shallow depth, so the effects of radiation penetration must be included for a realistic prediction of the ocean mixed layer. For this purpose, Paulson and Simpson’s (1977) formula type I was used, such as RSW = −I0[0.58 exp(−z/ζ1) + 0.42 exp(−z/ζ2), where I0 is the surface irradiance, ζ1 = 0.35 m, and ζ2 = 23 m.

Figure 9a shows the evolution of buoyancy profiles obtained from the simulation corresponding to EXP LW, where radiation penetration is included (EXP LW_R). As expected, a large amount of buoyancy is transferred downward below the mixed layer, which makes B(z = 0), corresponding to sea surface temperature (SST), substantially smaller than in Fig. 2. The thermocline also becomes slightly deeper because turbulence is less strongly suppressed by stratification. The general pattern of buoyancy profile evolution remains invariant, however. The time series of Rf, Ri, and r at z = 12.5 m follow a similar pattern, but their equilibrium values are slightly smaller (Fig. 10). It suggests that both N2 and S2 decrease at the thermocline in EXP LW_R, but Ri [= (N/S)2] tends to be smaller because radiation penetration allows the downward transport of buoyancy but not of momentum.

Fig. 9.

Evolution of vertical profiles of B with time (Δt = 1 h, and the dotted line represents t = 0.5 h) (a) in the presence of radiation penetration and (b) under the diurnally varying surface buoyancy flux.

Fig. 9.

Evolution of vertical profiles of B with time (Δt = 1 h, and the dotted line represents t = 0.5 h) (a) in the presence of radiation penetration and (b) under the diurnally varying surface buoyancy flux.

Fig. 10.

Time series of physical variables at z = 12.5 m (EXP LW, solid; EXP LW_R, dashed): (a) Rf, (b) Ri, and (c) /( + ).

Fig. 10.

Time series of physical variables at z = 12.5 m (EXP LW, solid; EXP LW_R, dashed): (a) Rf, (b) Ri, and (c) /( + ).

Also investigated was how the diurnal thermocline formation process is modified under the diurnal variation of the surface buoyancy flux. For this purpose the surface buoyancy flux was imposed in the form Q0[1−cos(2πt/T)] with T = 8 h, for the simulation corresponding to EXP LW. As expected, the initial thermocline formation is slower (t ∼ 2 h) because Q0 is smaller, but the stronger thermocline and higher B(z = 0) appear when Q0 is larger (t ∼ 4 h; Fig. 9b). It is interesting to observe that B(z = 0) decreases when Q0 decreases after t = 4 h in spite of the continuous supply of buoyancy into the mixed layer, because the depth of the thermocline increases. The corresponding decrease of SST is also observed from observation data shown in Fig. 3 (1646 LT).

g. Sensitivity to the vertical resolution

Unlike the previous LES dealing with the homogeneous mixed layer (Noh et al. 2004), stratification is generated in the present LES. The vertical motion of eddies is suppressed under stratification, and the results are more likely to be affected by the vertical resolution. Therefore, EXP LW was repeated with half the vertical grid size (Δz = 0.625 m; EXP LW_H).

The overall results shown in previous sections (Figs. 1, 2, 4, and 5) remain essentially invariant under EXP LW_H. Meanwhile, the vertical TKE at the thermocline is slightly larger in EXP LW_H, causing the slight decrease of Rf and Ri and the slight increase of r at the thermocline (Ri ≅ 0.25, Rf ≅ 0.3, and r ≅ 0.2; Fig. 11). It is also found that the transport of heat and momentum to below the mixed layer is slightly delayed, as can be inferred from the large fluctuation of Rf at z = 12.5 m during the initial stage in Fig. 11.

Fig. 11.

Time series of physical variables at z = 12.5 m (EXP LW, solid; EXP LW_H, dashed): (a) Rf, (b) Ri, and (c) /( + ).

Fig. 11.

Time series of physical variables at z = 12.5 m (EXP LW, solid; EXP LW_H, dashed): (a) Rf, (b) Ri, and (c) /( + ).

4. Conclusions

Large-eddy simulation of the ocean mixed layer under a stabilizing buoyancy flux was performed. The formation of a diurnal thermocline was reproduced successfully, in agreement with observation and mixed layer model results. The analysis of the simulation result revealed that the formation of a diurnal thermocline passes through two different phases: the formation of a thermocline (formation stage) and increasing thickness of the thermocline thereafter (growth stage). It was also found that the TKE flux divergence dominates TKE production within the mixed layer, but the local equilibrium PsPb − ɛ = 0 is achieved at the thermocline. The turbulence maintained by the local equilibrium at the thermocline causes stratification below the mixed layer. Based on these features, two different feedback mechanisms that control the formation and growth stages were suggested (Table 1).

Various physical variables at the thermocline—such as S2, N2, Rf, Ri, Km, Kh, vertical and horizontal TKE, and length scales—vary significantly during the formation stage, but they approach constant values in the growth stage thereafter. The development of S2 lags behind N2 by 0.5–1 h, causing the initial appearance of abnormally large values of Rf and Ri. The mixing length scale may be related to the buoyancy length scale lb(= q/N), but the dissipation length scale is not.

It was found that a diurnal thermocline cannot be formed in the absence of both WB and LC. In this case, TKE production is dominated by shear production near the sea surface, and strong stratification and shear appear near the surface, similar to the atmospheric boundary layer. Meanwhile, a diurnal thermocline is formed in the presence of either WB or LC only, although the depth of the thermocline tends to be shallower. Radiation penetration makes the depth of the thermocline slightly deeper and the values of Rf, Ri, and r at the thermocline slightly smaller.

The present work was focused on understanding the dynamical process of diurnal thermocline formation. For a realistic prediction of the response of the ocean mixed layer to the surface heating, however, it may be necessary to consider various other aspects, such as the appearance of strong stratification near the sea surface under calm weather and the generation, propagation, and breaking of surface gravity waves. It is expected that a similar mechanism underlies seasonal thermocline formation, but a separate study may be required for its understanding because of the much longer time scale and the diurnal cycle of heat flux. Meanwhile, further analysis of the LES data can be applied to elaborate the ocean mixed layer model, especially in the presence of WB and LC (e.g., Li and Garrett 1997; D’Alessio et al. 1998; Noh 2004; Kantha and Clayson 2004).

Acknowledgments

This work was funded by the Korea Meteorological Administration Research and Development Program under Grant CATER 2006-4201. Supercomputing was supported from KISTI.

REFERENCES

REFERENCES
Agrawal
,
Y. C.
,
E. A.
Terray
,
M. A.
Donelan
,
P. A.
Hwang
,
A. J.
Williams
III
,
W. M.
Drennan
,
K. K.
Kahma
, and
S. A.
Kitaigorodskii
,
1992
:
Enhanced dissipation of kinetic energy beneath surface waves.
Nature
,
359
,
219
220
.
Brainerd
,
K. E.
, and
M. C.
Gregg
,
1993
:
Diurnal restratification and turbulence in the oceanic surface mixed layer. 1. Observations.
J. Geophys. Res.
,
98
,
22645
22656
.
Britter
,
R. E.
,
J. C. R.
Hunt
,
G. L.
Marsh
, and
W. H.
Snyder
,
1983
:
The effects of stable stratification on turbulent diffusion and the decay of grid turbulence.
J. Fluid Mech.
,
127
,
27
44
.
Caldwell
,
D. R.
,
R-C.
Lien
,
J. N.
Moum
, and
M. C.
Gregg
,
1997
:
Turbulence decay and restratification in the equatorial ocean surface layer following nighttime convection.
J. Phys. Oceanogr.
,
27
,
1120
1132
.
Caughey
,
S. J.
,
J. C.
Wyngaard
, and
J. C.
Kaimal
,
1979
:
Turbulence in the evolving stable boundary layer.
J. Atmos. Sci.
,
36
,
1041
1052
.
Craig
,
P. D.
, and
M. L.
Banner
,
1994
:
Modeling wave-enhanced turbulence in the ocean surface layer.
J. Phys. Oceanogr.
,
24
,
2546
2559
.
Craik
,
A. D. D.
, and
S.
Leibovich
,
1976
:
A rational model for Langmuir circulations.
J. Fluid Mech.
,
73
,
401
426
.
D’Alessio
,
J. S. D.
,
K.
Abdella
, and
N. A.
McFarlane
,
1998
:
A new second-order turbulence closure scheme for modeling the ocean mixed layer.
J. Phys. Oceanogr.
,
28
,
1624
1641
.
Davis
,
R. E.
,
R.
de Szoeke
, and
P. P.
Niiler
,
1981
:
Variability in the upper ocean during MILE. Part II: Modeling the mixed layer response.
Deep-Sea Res.
,
28A
,
1453
1475
.
Delnore
,
V. E.
,
1972
:
Diurnal variations of temperature and energy budget for the ocean mixed layer during BOMEX.
J. Phys. Oceanogr.
,
2
,
239
247
.
Drennan
,
W. M.
,
M. A.
Donelan
,
E. A.
Terray
, and
K. B.
Katsaros
,
1996
:
Oceanic turbulence dissipation measurements in SWADE.
J. Phys. Oceanogr.
,
26
,
808
815
.
Fairall
,
C. W.
,
E. F.
Bradley
,
J. E.
Hare
,
A. A.
Grachev
, and
J. B.
Edson
,
2003
:
Bulk parameterization of air–sea fluxes: Updates and verification for the COARE algorithm.
J. Climate
,
16
,
571
591
.
Gemmrich
,
J. R.
, and
D. M.
Farmer
,
1999
:
Near-surface turbulence and thermal structure in a wind-driven sea.
J. Phys. Oceanogr.
,
29
,
480
499
.
Hopfinger
,
E. J.
, and
P. F.
Linden
,
1982
:
Formation of thermoclines in zero-mean-shear turbulence subjected to a stabilizing buoyancy flux.
J. Fluid Mech.
,
114
,
157
173
.
Hunt
,
J. C. R.
,
D. D.
Stretch
, and
R. E.
Britter
,
1988
:
Length scales in stably stratified turbulent flows and their use in turbulence models.
Stably Stratified Flow and Dense Gas Dispersion, J. S. Puttock, Ed., Oxford University Press, 285–322
.
Kantha
,
L. H.
, and
R. R.
Long
,
1980
:
Turbulent mixing with stabilizing surface buoyancy flux.
Phys. Fluids
,
23
,
2142
2143
.
Kantha
,
L. H.
, and
C. A.
Clayson
,
2000
:
Small-Scale Processes in Geophysical Fluid Flows.
Academic Press, 888 pp
.
Kantha
,
L. H.
, and
C. A.
Clayson
,
2004
:
On the effect of surface gravity waves on mixing in the ocean mixed layer.
Ocean Modell.
,
6
,
101
124
.
Kraus
,
E. B.
,
1988
:
Merits and demerits of different approaches to mixed layer modelling.
Small-Scale Turbulence and Mixing in the Ocean, J. C. J. Nihoul and B. M. Jamart, Eds., Elsevier, 37–50
.
Kudryavtsev
,
N. V.
, and
A. V.
Soloviev
,
1990
:
Slippery near-surface layer of the ocean arising due to daytime solar heating.
J. Phys. Oceanogr.
,
20
,
617
628
.
Li
,
M.
, and
C.
Garrett
,
1997
:
Mixed layer deepening due to Langmuir circulation.
J. Phys. Oceanogr.
,
27
,
121
132
.
Li
,
M.
,
C.
Garrett
, and
E.
Skyllingstad
,
2005
:
A regime diagram for classifying large eddies in the upper ocean.
Deep-Sea Res. I
,
52
,
259
278
.
McWilliams
,
J. C.
,
P. P.
Sullivan
, and
C. H.
Moeng
,
1997
:
Langmuir turbulence in the ocean.
J. Fluid Mech.
,
334
,
1
30
.
Min
,
H. S.
, and
Y.
Noh
,
2004
:
Influence of the surface heating on Langmuir circulation.
J. Phys. Oceanogr.
,
34
,
2630
2641
.
Niiler
,
P. P.
, and
E. B.
Kraus
,
1977
:
One-dimensional models of the upper ocean.
Modelling and Prediction of the Upper Layers of the Ocean, E. B. Kraus, Ed., Pergamon, 143–172
.
Noh
,
Y.
,
1996
:
Dynamics of diurnal thermocline formation in the oceanic mixed layer.
J. Phys. Oceanogr.
,
26
,
2183
2195
.
Noh
,
Y.
,
2004
:
Sensitivity to wave breaking and the Prandtl number in the ocean mixed layer model and its dependence on latitude.
Geophys. Res. Lett.
,
31
,
L23305
.
doi:10.1029/2004GL021289
.
Noh
,
Y.
, and
R. R.
Long
,
1990
:
Turbulent mixing in a rotating, stratified fluid.
Geophys. Astrophys. Fluid Dyn.
,
53
,
125
143
.
Noh
,
Y.
, and
H. J. S.
Fernando
,
1991
:
A numerical study on the formation of a thermocline in shear-free turbulence.
Phys. Fluids
,
3A
,
422
426
.
Noh
,
Y.
, and
H. J.
Kim
,
1999
:
Simulations of temperature and turbulence structure of the oceanic boundary layer with the improved near-surface process.
J. Geophys. Res.
,
104
,
15621
15634
.
Noh
,
Y.
,
H. S.
Min
, and
S.
Raasch
,
2004
:
Large-eddy simulation of the ocean mixed layer: The effects of wave breaking and Langmuir circulation.
J. Phys. Oceanogr.
,
34
,
720
735
.
Noh
,
Y.
,
I. S.
Kang
,
M.
Herold
, and
S.
Raasch
,
2006
:
Large-eddy simulation of particle settling in the ocean mixed layer.
Phys. Fluids
,
18
,
085109
.
doi:10.1063/1.2337098
.
Paulson
,
C. A.
, and
J. J.
Simpson
,
1977
:
Irradiance measurements in the upper ocean.
J. Phys. Oceanogr.
,
7
,
952
956
.
Polton
,
J. A.
, and
S. E.
Belcher
,
2007
:
Langmuir turbulence and deeply penetrating jets in an unstratified mixed layer.
J. Geophys. Res.
,
112
,
C09020
.
doi:10.1029/2007JC004205
.
Raasch
,
S.
, and
M.
Schröter
,
2001
:
PALM—A large-eddy simulation model performing on massively parallel computers.
Meteor. Z.
,
10
,
363
372
.
Schumann
,
U.
, and
T.
Gerz
,
1995
:
Turbulent mixing in stably stratified shear flows.
J. Appl. Meteor.
,
34
,
33
48
.
Skyllingstad
,
E. D.
, and
D. W.
Denbo
,
1995
:
An ocean large-eddy simulation of Langmuir circulations and convection in the surface mixed layer.
J. Geophys. Res.
,
100
,
8501
8522
.
Skyllingstad
,
E. D.
,
W. D.
Smyth
,
J. N.
Moum
, and
H.
Wijesekera
,
1999
:
Upper-ocean turbulence during a westerly wind burst: A comparison of large-eddy simulation results and microstructure measurements.
J. Phys. Oceanogr.
,
29
,
5
28
.
Soloviev
,
A.
, and
R.
Lukas
,
1997
:
Observation of large diurnal warming events in the near-surface layer of the western equatorial Pacific warm pool.
Deep-Sea Res. I
,
44
,
1055
1076
.
Stommel
,
H.
,
K.
Saunders
,
W.
Simmons
, and
J.
Cooper
,
1969
:
Observation of the diurnal thermocline.
Deep-Sea Res.
,
16
,
269
284
.
Sullivan
,
P. P.
,
J. C.
McWilliams
, and
W. K.
Melville
,
2007
:
Surface gravity wave effects in the oceanic boundary layer: Large-eddy simulation with vortex force and stochastic breakers.
J. Fluid Mech.
,
593
,
405
452
.
Woods
,
J. D.
, and
W.
Barkmann
,
1986
:
The response of the upper ocean to solar heating. I. The mixed layer.
Quart. J. Roy. Meteor. Soc.
,
112
,
1
27
.

Footnotes

Corresponding author address: Yign Noh, Department of Atmospheric Sciences/Global Environmental Laboratory, Yonsei University, 134 Shinchon-dong, Seodaemun-gu, Seoul, South Korea. Email: noh@yonsei.ac.kr