Large-eddy simulation (LES) is used to investigate how dominant breaking waves in the ocean under hurricane-force winds affect the drag and near-surface airflow turbulence. The LES explicitly resolves the wake turbulence produced by dominant-scale breakers. Effects of unresolved roughness such as short breakers, nonbreaking waves, and sea foam are modeled as the subgrid-scale drag. Compared to the laboratory conditions previously studied using the same method, dominant-scale breakers in open-ocean conditions are less frequent, and the subgrid-scale drag is more significant. Nevertheless, dominant-scale breakers are more fully exposed to high winds and produce more intense wakes individually. As a result, they support a large portion of the total drag and significantly influence the turbulence for many ocean conditions that are likely to occur. The intense wake turbulence is characterized by flow separation, upward bursts of wind, and upward flux of the turbulent kinetic energy (TKE), all of which may influence sea spray dispersion. Similarly to the findings in the laboratory conditions, high production of wake turbulence shortcuts the inertial energy cascade, causes high TKE dissipation, and contributes to the reduction of the drag coefficient. The results also indicate that if the drag coefficient decreases with increasing wind at very high winds, as some recent observations suggest, then the unresolved roughness must also decrease.
Air–sea exchanges of momentum, heat, and gas as well as the suspension and dispersion of water droplets and other passive tracers are key factors affecting many atmospheric and oceanic processes. These factors are regulated by airflow turbulence near the air–sea boundary. Despite its importance, such turbulence is poorly understood at high winds, because it is affected by complex physical processes such as sea foam (Powell et al. 2003; Soloviev and Lukas 2010; Holthuijsen et al. 2012), sea spray (Makin 2005; Barenblatt et al. 2005; Bianco et al. 2011; Bao et al. 2011; Kudryavtsev and Makin 2011), and breaking waves (Kudryavtsev and Makin 2007; Kukulka et al. 2007; Kukulka and Hara 2008b; Mueller and Veron 2009; Banner and Morison 2010). Breaking waves at high winds induce wake turbulence, which is distinctively different from the regular shear turbulence (Reul et al. 2008; Suzuki et al. 2013). One of the wake turbulence features is airflow separation (Reul et al. 2008; Suzuki et al. 2013). The airflow separation over a breaker affects the form drag over the breaker itself (Kudryavtsev and Makin 2007; Kukulka et al. 2007; Kukulka and Hara 2008b; Mueller and Veron 2009). Moreover, the separated flow region shelters smaller-scale waves from direct wind forcing (sheltering effect) and further modifies the air–sea momentum flux (Kudryavtsev and Makin 2007; Kukulka et al. 2007; Kukulka and Hara 2008b; Mueller and Veron 2009).
The effects of breakers and wake turbulence at high winds have been demonstrated using large-eddy simulation (LES) by Suzuki et al. (2013) for laboratory-scale breakers that are short (i.e., wavelength <1 m) and narrowbanded. In their study, the impact of intermittent and transient wave breaking events is modeled as localized form drag, which generates wake turbulence including airflow separation. They find that more than 90% of the total air–sea momentum flux is due to the form drag of breakers; that is, the contributions of the nonbreaking wave form drag and the surface viscous stress are small. A similar result was obtained for laboratory conditions by Kudryavtsev and Makin (2007). Suzuki et al. (2013) also find that breaker form drag impedes the shear production of the turbulent kinetic energy (TKE) near the surface and, instead, produces a large amount of small-scale wake turbulence by transferring energy from large-scale motions (such as mean wind and gusts). This process shortcuts the inertial energy cascade and results in large TKE dissipation (integrated over the surface layer) and a smaller drag coefficient.
Although their study highlights the importance of breakers, their results are limited to laboratory conditions. In the open ocean, waves are longer and broadbanded, and the breaking of dominant waves near the spectral peak (hereafter, dominant-scale breaker) is relatively rare (Banner et al. 2000, 2002; Gemmrich et al. 2008; Kleiss and Melville 2010). It is often suggested that waves much shorter than the dominant waves support most of the total drag at the ocean surface and that the contribution of dominant-scale breakers is not significant at low-to-moderate wind speeds (Mueller and Veron 2009; Kudryavtsev and Makin 2007). However, an individual dominant-scale breaker may support a large form drag (Babanin et al. 2007), and it is unclear to what degree dominant-scale breakers affect the air–sea momentum flux and airflow turbulence at very high wind speeds.
To address this problem, we extend the LES approach of Suzuki et al. (2013) to open-ocean conditions in high winds. The strength of this approach is that the wake turbulence generated by dominant-scale breakers is explicitly simulated. This is in contrast to other model studies (viz., Kudryavtsev and Makin 2007; Kukulka et al. 2007; Kukulka and Hara 2008b; Mueller and Veron 2009) where the effects of wakes are simply parameterized. In this study we do not attempt to predict the drag coefficient, since the drag coefficient strongly depends on processes unresolved by the LES. Instead, the main purpose of this study is to understand to what degree dominant-scale breakers may affect the drag and airflow turbulence in various ocean conditions.
The ocean surface under hurricane-force winds has breaking and nonbreaking waves over a wide range of scales as well as other small-scale roughness, such as sea foam. An ideal numerical method of studying the sea surface influences on the airflow would be to simulate both the airflow and the sea surface, resolving all relevant scales. However, because a direct numerical simulation involving a high Reynolds number and broadband breaking and nonbreaking waves in a large domain is not a viable option [for a review, see Perlin et al. (2013)], we will not simulate the waves directly. Instead, we adopt a simpler LES approach with modeled wave effects.
LES is a numerical technique that has a limited resolution. It explicitly simulates motions at the resolved scales, but physical processes at the unresolved scales must be parameterized. In particular, because of the finite vertical resolution, our LES cannot resolve impacts of waves below a certain scale. The wave-induced stress (Belcher 1999; Kudryavtsev and Makin 2001; Kukulka and Hara 2008a; Mueller and Veron 2009; Banner and Morison 2010)—or the pressure flux in a surface-fitted coordinate system (Sullivan et al. 2000)—modifies the airflow turbulence above the water surface. Although there is some uncertainty, the height over which this wave influence is significant, roughly scales with the inner-layer height for nonbreaking waves (Belcher 1999) and with the breaker amplitude for breaking waves (Kukulka et al. 2007; Kukulka and Hara 2008a,b; Mueller and Veron 2009). Because the scale of ocean waves spans a wide range, the influence of waves below a certain scale falls below the LES vertical resolution and must be treated as part of the subgrid-scale (SGS) bottom boundary stress. In contrast, the effect of larger waves may be modeled as momentum injection within the LES domain and be explicitly resolved.
The vertical resolution of this study is chosen in such a way that it allows explicit simulations of wake turbulence made by dominant-scale breakers and, yet, is coarse enough to keep the computational cost feasible. To define a suitable vertical resolution, we make the following two assumptions according to Kukulka et al. (2007), Kukulka and Hara (2008a,b), and Mueller and Veron (2009):
The inner-layer height and the breaker amplitude are related to the wavenumber as δ/k and s/k, respectively, where k is the wavenumber and the nondimensional parameters δ and s range as 0.01 ≲ δ ≲ 0.1 (Kukulka and Hara 2008a,b) and 0.15 ≲ s ≲ 0.55 (Perlin et al. 2013). In particular, we assume that δ is roughly 0.05 and s = 0.3, following Kukulka and Hara (2008a,b).
The drag supported by nonbreaking waves longer than the spectral peak wavelength is negligible; hence, the longest nonbreaking waves relevant to our problem are those at the spectral peak.
With these assumptions, we may set the LES vertical resolution to be coarse enough to treat the form drag of all relevant nonbreaking waves as part of the subgrid-scale bottom boundary stress. At the same time, the vertical resolution is set fine enough to treat the form drag of dominant-scale breakers as momentum injections at resolved heights. In this way, we can explicitly simulate the airflow disturbances directly generated by dominant-scale breakers. Simulations of these disturbances and the form drag of dominant-scale breakers are detailed in sections 2b and 2c.
The subgrid-scale bottom boundary stress represents the net drag due to all relevant nonbreaking waves, unresolved-scale breakers, and other unresolved-scale roughness such as sea foam. In principle, it is possible to model the subgrid-scale drag in a sophisticated way by considering particular conditions of the wave spectrum, short breakers, and sea foam. However, in this study, we do not use such a model. Instead, the combined effect is represented using a single parameter called the drag coefficient of the unresolved roughness. This is because our knowledge of the wave spectrum, short breakers, and sea foam is very limited at hurricane-force winds. This simple approach also allows us to consider variations of these roughness elements as they naturally vary due to variable environmental conditions (e.g., wind history, ocean currents, fetch, and surfactant). We assume that the resolved-scale dynamics depends only on the net effect of the unresolved roughness elements, not on the details of wave spectrum or sea foam. The parameterization of the subgrid-scale bottom boundary stress is detailed in section 2d.
b. LES model of the boundary layer with breaker effects
In this study, only the atmospheric side of the surface layer is simulated. To include the breaker effect, our airflow simulations are coupled with simulations of breaker positions that evolve with time. The LES equations used are identical to those of Suzuki et al. (2013). The only difference between Suzuki et al. (2013) and the current study is the surface wave fields (sea states) considered. Thus, we review our LES framework only briefly.
When the wind blows over and around a breaker, a localized pressure perturbation appears at the air–sea interface and in the interior of the air surrounding the breaker. This pressure perturbation causes exchanges of momentum and energy between the breaker and the surrounding air. It is these momentum and energy exchanges that are modeled in our LES as the effects of breakers. The effects of the surface undulation and surface orbital velocities of breakers are assumed to be secondary. Thus, the bottom boundary (i.e., the air–sea interface) is idealized as a flat surface, and the LES equations are approximated with their Cartesian forms. To model the momentum and energy exchanges without the actual surface undulation, forcing terms are added in the LES equations. Although such a representation of the breaker effect is highly idealized, a flow simulation over and around a breaker using this modeling reproduces realistic flow separation and other essential characteristics of the wake turbulence (Suzuki et al. 2013).
where filtered variables are denoted by overbars; x1, x2, and x3 (or equivalently x, y, and z) are the streamwise, spanwise, and vertical coordinates, respectively; (u1, u2, u3) = (u, υ, w) are the velocity components in (x1, x2, x3) = (x, y, z), respectively; p is the pressure divided by a constant air density; is an external large-scale forcing used to drive the flow ( is constant in time, uniform in space, and positive); is the SGS stress; is the SGS kinetic energy; TSGS is the SGS transport; ϵ is the viscous dissipation; and t is time. The effects of a local discrete breaking wave event m are represented by and Wm, where the former is the momentum input to the resolved motion and the latter is work done to the SGS turbulence. As our focus is on the effects of breakers in a relatively thin roughness sublayer, the effect of stratification and the Coriolis force are not considered. Although the effects of sea sprays and sea foam are not explicitly included in the LES equation, their effects are implicitly included in the surface stress boundary condition as discussed below.
In Eqs. (1) and (3), the regular SGS terms (viz., Rij, TSGS, and ϵ) and the breaker effect terms (viz., and Wm) require modeling. The regular SGS terms are modeled using a conventional TKE closure SGS parameterization describe by Moeng (1984). In Suzuki et al. (2013), another TKE closure SGS parameterization described by Sullivan et al. (1994) was also tested, and they found only minor differences in the results between the two parameterizations. The SGS stress is modeled with eddy viscosity νT diagnosed based on e; TSGS is modeled as the downgradient diffusion of e, namely, (∂/∂xj)(2νT∂e/∂xj); and ϵ is assumed to be proportional to e3/2.
The forcing due to the breaker-induced pressure perturbation on the airflow around a breaker m is modeled based on a conventional aerodynamic form drag formula (e.g., Kukulka et al. 2007). The forcing is localized and appears only inside a volume surrounding the breaker m. The dimensions of the volume are empirically determined to be height × along-crest length × across-crest length = a × λ × 0.5λ, where λ and a are the wavelength and amplitude of the breaker, respectively. Within this volume,
where H = 2a is the height of the frontal area of the breaker; is an empirically determined form drag coefficient of the breaker and is assumed to be the same for all breakers; c is the propagation velocity of the breaker and is assumed to be related to the (angular) wavenumber k and the gravitational acceleration g by ; and uAT is set to the instantaneous wind parallel to c at the upwind side of the breaker and at a height z = a. If uAT is opposite to or slower than c, then is set to zero. As mentioned in section 2a, we assume that ak = s where s = 0.3. Outside the volume surrounding the breaker m, the momentum input .
The work done on the SGS turbulence by the pressure perturbation around the breaker m is modeled as
for the following reason. First, we may derive the equation for the resolved energy by multiplying Eq. (1) and . In this equation, the rate of work done by breaker forcing on resolved winds is given by . On the other hand, the rate of energy transfer to the breaker may be estimated by the breaker propagation velocity times the form drag, namely, . Then, the conservation of energy may be written as . To satisfy the energy conservation, we simply model Wm as . When Wm is negative, we reset at that location to avoid an unphysical SGS work input (i.e., breakers do not convert SGS motions into resolved-scale motions).
c. Field of breakers
During our LES runs, positions of discrete breakers over a range of wavenumbers are generated intermittently in time, randomly in space, and independently from the airflow. Once generated, each breaker lasts for one wave period , and its position moves at its breaker propagation velocity c. The crest length of each breaker is set to its wavelength. These parameter choices follow Sullivan et al. (2007) and Suzuki et al. (2011, 2013), and our results are relatively insensitive to the particular choices made here. In our simulations, a random number of breakers at each wavenumber are generated at each time step in such a way that the resultant breaker field satisfies a specified breaking wave distribution Λ.
The breaking wave distribution Λ is defined such that Λ(k, σ)kdkdσ represents the average length of breaking crests per unit horizontal area of the sea surface for waves with their wavenumbers between k − dk/2 and k + dk/2 and their propagation directions between σ − dσ/2 and σ + dσ/2 (e.g., Phillips 1985; Kleiss and Melville 2011). It is also common to use the one-dimensional breaking wave distribution defined as
The distribution Λ(k) may be converted from or to Λ(c) based on Λ(k)dk/dc = Λ(c) and an empirical relationship like (Melville and Matusov 2002). In this study, however, we use a linear wave relationship for simplicity. Note that Λ(k) is dimensionless whereas Λ(c) is not.
Our breaking wave distributions are specified according to observational data at the open ocean. Unfortunately, direct observations of Λ are available for winds only up to U10 ≈ 20 m s−1 (Kleiss and Melville 2010), where U10 is the mean wind speed at 10-m height. Thus, we estimate Λ at higher winds with the aid of observational whitecap coverage data, which are available for winds up to U10 ≈ 50 m s−1 (Holthuijsen et al. 2012). Since whitecap coverage is related to Λ and may be proportional to (e.g., Kleiss and Melville 2010), it can be used to deduce approximate Λ at U10 > 20 m s−1. First, let WDB be the whitecap coverage due to the dominant-scale breakers, whose wavenumbers range between kd and nkd (or the phase speed between cd and ), where n is a number greater than 1. We assume that the breakers longer than λd = 2π/kd occur so rarely that they have negligible contributions to the whitecap coverage. Then, WDB may be expressed as
where λ′ = kdλ and k′ = k/kd are the normalized wavelength and normalized wavenumber. The rhs of Eq. (7) is expressed all in dimensionless quantities. Hence, Eq. (7) implies that WDB changes with change in Λ(k′) but not with change in the dominant scale kd.
A recent observation by Holthuijsen et al. (2012) shows no systematic dependency of whitecap coverage on the wind speed for U10 > 24 m s−1. Their measurements of the whitecap coverage fluctuate around 4% and are always below 10%. Note that whitecap coverage at U10 ≈ 20 m s−1 is also generally between a few percent and 10% at most (e.g., Kleiss and Melville 2010; Holthuijsen et al. 2012). We assume that the observed whitecap coverage is mainly related to the distribution of dominant-scale breakers and, hence, is approximately WDB. Then, to keep WDB independent of the wind speed, the overall value of Λ(k′) in Eq. (7) has to be independent of the wind speed; that is, the overall value of Λ(k′) at high winds should be similar to that at U10 ≈ 20 m s−1, where more detailed data of Λ are available. [Note that constancy of WDB implies constancy of the overall value of Λ(k′), but does not imply constancy of Λ(c). In fact, the value of Λ(c) strongly depends on kd even if WDB is constant, as will be shown in Fig. 1.] Note that the functional shape of Λ is not as important as the overall value of Λ (Suzuki et al. 2013). Therefore, in this study, we consider a Λ(k′) directly measured at U10 ≈ 20 m s−1 as the baseline Λ(k′) and use the same Λ(k′) at all higher wind speeds. Then, we consider some variations of Λ around this baseline. We also consider variation in the dominant scale of the breakers: namely, kd or equivalently λd. Hereafter, the baseline Λ is denoted by ΛB.
In particular, to obtain ΛB, we use a Λ observed at U10 = 17 m s−1 by Kleiss and Melville (2010) (Fig. 1). When this Λ was observed, the spectral peak wavelength λp was 29.22 m, and the wave age was 9, where cp is the phase speed at the spectral peak and is the friction velocity. We set λd in this observation equal to the observed λp. Assigning a longer λd does not change our results since breakers longer than this λp occur very infrequently, according to the reported Λ.
There are many other breaking wave distributions shown in Kleiss and Melville (2010). However, the functional shapes of the other distributions are similar to this one in the sense that they can be roughly approximated by multiplying a constant to this particular Λ and setting an appropriate dominant scale. Based on this trend, we assume that it is sufficient to specify various distributions Λ(k) by simply multiplying a variable constant β to ΛB(k/kd) and varying kd.
In our LES, there is a cutoff wavenumber ko to truncate a specified Λ(k), where ko is the wavenumber of the largest resolved breaker considered. The corresponding wavelength is λo = 2π/ko. After testing, we set ko = kd or equivalently λo = λd. Assigning a longer λo does not change our results since breakers longer than λd are very rare in the conditions considered in this study. For the shortest resolved breakers, Λ is truncated at 9ko or equivalently λo/9. This is because the breaker amplitude at this cutoff wavenumber corresponds to the height of our first LES vertical grid level. Breakers shorter than this scale are treated as unresolved-scale roughness.
In addition to Λ(k), our breaker field needs a specification of the directional distribution of the breakers. The directional distribution of Λ in Kleiss and Melville (2010) is roughly cos3ασ where 0.5 < α < 3, depending on the measurement methods. In our study, however, we focus on a special case where the dominant-scale breakers are unidirectional; that is, all breaking crests of the dominant-scale breakers are perpendicular to the mean wind. This is because, according to our preliminary simulations, the impact of a breaking wave that propagates at an angle σ relative to the mean wind direction is practically (on average) identical to that of the same breaker that propagates in the mean wind direction but with the form drag coefficient in Eq. (4) reduced by a factor . Therefore, introducing a wider directional distribution is practically equivalent to slightly reducing in the unidirectional case. As discussed below, we examine a wide range of in this study.
d. Unresolved roughness parameterization
In the previous subsection, we modeled the field of dominant-scale breakers, whose individual occurrences are resolved in our simulations. In this subsection, we model all the other types of roughness, which are not explicitly resolved. Such roughness includes unresolved breakers, all nonbreaking waves, viscous stress, and other processes such as sea foam and sea spray.
where i = 1, 2; κ is the von Kármán constant; z1 is the height of the lowest grid level for ; is the roughness length for the unresolved roughness; and U1 is the mean wind speed at z = z1. At the bottom boundary, the SGS stress Ri3 is replaced with . This parameterization assumes that the unresolved roughness is isotropic; that is, the bottom drag is always aligned with and opposite to the local wind. However, the roughness of the actual ocean surface is likely to be anisotropic as the waves and the form drag over them have some particular directionality. For example, at a very young wave age where the directional spreading is narrow, the drag may be highly anisotropic and be idealized as , while (i.e., the stress in the mean wind direction) is still given by Eq. (8) (Suzuki et al. 2011). At the open ocean, the directionality of the unresolved roughness lies somewhere between the isotropic and highly anisotropic conditions. In the current study, we simply use the highly anisotropic parameterization for the following reasons:
When the isotropic formula is used instead, the results of the drag coefficient CD10 are slightly reduced (Suzuki et al. 2011), where . However, such difference is systematic, and our conclusions essentially remain the same.
Unlike the isotropic parameterization, the highly anisotropic parameterization suffers little from a known deficiency in LES of boundary layers (Suzuki et al. 2011).
The unresolved roughness length can be expressed in terms of defined as
This expression is more convenient because is equal to CD10 when there are no resolved-scale breakers.
e. Sea state parameters and their ranges
In this subsection, we summarize the model input parameters and discuss their ranges. As previously discussed, we specify Λ(k/ko) by multiplying a constant factor β to the baseline ΛB(k/ko). Therefore, the sea state in our LES is parameterized with four parameters: namely, β, λo = 2π/ko, , and . The first parameter β controls the overall amount of the dominant-scale breakers. The second parameter λo is the wavelength of the largest resolved breakers. It also represents the dominant scale of the breakers as ko = kd. The third parameter is the aerodynamic form drag coefficient for individual breakers and relates the local wind forcing on a breaker and the resulting form drag. The fourth parameter is the bulk drag coefficient for unresolved roughness.
After preliminary testing, we have found, at all conditions considered in this study, that the effect of increasing by a certain factor is practically the same as the effect of increasing β by the same factor. Thus, by considering the product quantity , we may reduce the number of the sea state parameters to three: namely, (hereafter called the breaker factor), λo, and .
To constrain the values of these input parameters, we rely on previous observational and theoretical studies. There is a general consensus that the value of is of order one (Kukulka et al. 2007; Kudryavtsev and Makin 2001; Mueller and Veron 2009). In addition, Suzuki et al. (2013) estimated the range of to be roughly between 0.5 and 3.0 by comparing their large-eddy simulations to an observation (Reul et al. 2008) of breaker-induced wakes. An open ocean undergoes various combinations of wind and sea states. The breaking distribution Λ is known to vary, at least in moderate winds, by an order of magnitude (e.g., Kleiss and Melville 2010) at a given wind speed. Thus, a quite wide variation in the breaker factor may be expected. However, we do not expect Λ of the dominant-scale breakers to reduce with increasing wind; that is, β is not expected to be much less than 1 at high wind speeds. Thus, we do not explore fb less than 0.5. Our results also indicate that fb = 6 results in sufficiently large CD10 compared to the observational upper bound of CD10 at high wind speeds. Therefore, the range of fb explored in this study is between 0.5 and 6.
The range of explored is between 0 and 0.0022. As shown later, this covers most of the conditions that are expected to occur at the open ocean. For the lower bound of , we consider the possibility of the unresolved roughness becoming as smooth as a flat surface or even a slip surface. Flat walls have , where is the roughness length normalized with the wall friction velocity and the kinematic viscosity. For the wall friction velocity resulted in this study, is equivalent to , according to Eq. (9). At very high wind speeds, there are suggestions that a slip sea surface made of water–air emulsion may develop due to sea foam, bubbles, and spray (Powell et al. 2003; Soloviev and Lukas 2010; Holthuijsen et al. 2012). Therefore, we allow to decrease to zero. At U10 = 17 m s−1, we do not consider less than 0.0006 as a slip sea surface would not develop at such a wind condition.
The dominant scale of breakers is also variable at the open ocean. As discussed earlier, at U10 = 17 m s−1 we set λo = 29.22 m based on the observation. At the highest wind speed of U10 = 59 m s−1 examined in this study, rough estimates of the largest spectral peak wavelength and significant wave height under typical hurricane conditions are 291 and 14 m, respectively [according to Eqs. (10), (11), and (12) of Young (2003)]. We therefore consider two different dominant scales (lower and upper bounds): namely, λo = 29.22 m and λo = 256.46 m. The corresponding phase speeds for deep-water waves are 6.75 and 20.00 m s−1, respectively. For both cases, the wavelength of the shortest resolved breakers is λo/9; namely, 3.25 m for the shorter dominant scale and 28.50 m for the longer dominant scale. The corresponding range of fbΛB is shown in Fig. 1 (thin solid and dashed lines).
f. Numerical method
Time integration uses an explicit, third-order, three-substep Runge–Kutta scheme. A fixed time step is used based on a fixed Courant–Fredrichs–Lewy condition; = 0.063 to 0.26, depending on the simulated cases where ao = s/ko is the amplitude of the largest resolved breaker and the breaking threshold slope s = 0.3. Horizontal differentiation uses the pseudospectral method. Vertical differentiation uses the second-order-centered finite difference method on a vertically staggered grid. The variables , e, and Wm are stored at the same grid levels (hereafter, w-nodes), and , , , , and are stored at the grid levels (hereafter, u-nodes) located midway between the w-nodes. The w-nodes hold the bottom and top boundaries. The bottom boundary is at z = 0. The grid spacing is horizontally uniform and vertically nonuniform. We locate the fifth u-node at the height of the tallest breaker’s amplitude ao and set the distances of the lowest six w spacing to be Δz/ao = . Above this, each w-spacing Δz/ao is 1.03 times larger than the spacing one-node below. The horizontal boundaries are periodic. The top boundary is frictionless and nonpermeable. The bottom boundary is nonpermeable. The horizontal domain size is Lx × Ly, where Lx/ao = Ly/ao = 83.78; that is, Lx is 4 times λo. The domain height Lz is set such that Lz/ao = 56.22. The grid has 128 × 128 nodes horizontally and 96 nodes vertically.
The initial condition is a small and uniform streamwise wind everywhere. In reality, the wave field evolves in time or space; as a result, the airflow turbulence in such conditions may not be horizontally homogeneous or steady. However, in this study, we assume that the sea state (fb, , and λo) is constant in time and space (i.e., the wave growth in time or space is ignored). All data are obtained after the airflow has converged to a statistically steady (i.e., fully developed) state.
Some quantities are averaged for the following analysis. The averaging is done over each horizontal plane and over a long time to make sure that the errors in the momentum budget and TKE budget are less than 5% (viz., 18–81 large-eddy turnover time and 150–1260 wave periods of the longest breakers, depending on the simulation. Note that the averaging time can be reduced by increasing the horizontal domain size).
In the following, angle brackets denote a horizontal average, and a single prime denotes the deviation from it; for example, .
a. Drag coefficient CD10
In all simulations, the roughness sublayer appears adjacent to the modeled ocean surface (i.e., bottom surface), and the log layer establishes above the roughness sublayer. The results of the roughness length and the corresponding CD10 are computed by extending the log profile to the surface. We do not simply use the wind speed at 10-m height to calculate the drag coefficient because 10-m height may be inside the roughness sublayer where the wind profile is not logarithmic.
As expected, CD10 obtained depends on the wind speed U10 and the sea state is described by the breaker factor fb, the bulk drag coefficient for unresolved roughness , and the dominant scale of the breakers λo. Figure 2 presents the result of CD10 in several ways to show the dependence of CD10 on both U10 and the sea state parameters (fb, , and λo): that is, CD10 as a function of U10 (Fig. 2a), and CD10 as a function of fb and for three different combinations of U10 and λo (Figs. 2b,c,d). These figures also show the typical ranges of CD10 observed in the open ocean. At U10 ≈ 17 m s−1, the observational data of CD10 vary from 0.0010 ± 0.0001 (Holthuijsen et al. 2012), which is the lowest reported, to 0.0022 ± 0.0003 (Petersen and Renfrew 2009), which is the highest reported. At U10 ≈ 40 m s−1, the observational data vary from 0.0012 ± 0.0003 to 0.0030 ± 0.0025 (Vickery et al. 2009), and at U10 ≈ 59 m s−1 they vary from 0.0007 ± 0.0003 (Holthuijsen et al. 2012) to 0.0036 ± 0.0015 (Bell et al. 2012). In addition to CD10, Figs. 2b, 2c, and 2d show the percentage of the drag supported by the dominant-scale (resolved scale) breakers. Note, as discussed in section 2e, that corresponds to a condition that the unresolved roughness is equivalent to a flat surface. At U10 = 59 m s−1, the simulations have been extended to a slip surface; that is, .
The range of sea states explored in this study results in CD10 variations roughly comparable to the CD10 variations in the observational data. Figure 2a shows the values of CD10 at some of the roughest sea states and smoothest sea states considered in this study. These values of CD10 fall relatively close to the upper and lower values of the observational data. For reference, the sea states (i.e., the combinations of fb, , and λo) yielding these values of CD10 are indicated in Figs. 2b, 2c, and 2d (small and large circles and squares).
The wind dependence shown in Fig. 2a indicates that, when the sea state parameters are smaller and kept unchanged as the wind speed increases, CD10 increases to an upper-bound value and remains constant or saturates thereafter. It never decreases with further wind increase. A similar conclusion holds true for very short breakers in the laboratory, where breaking is very frequent (Donelan et al. 2004; Suzuki et al. 2013). The values of CD10 at the roughest sea state (fb, , and λo) = (6, 0.0022, and 256.46 m) do not show saturation within the wind speed considered (large circle in Fig. 2a).
When the wind is fixed, CD10 in high wind monotonically increases as fb or increases (Figs. 2b,c,d). The effect of increasing λo is also to increase CD10. However, this effect is negligibly small when fb and are small.
Our results of CD10 may be qualitatively consistent with a recent observation of the drag coefficient and sea states at high winds. Holthuijsen et al. (2012) observed CD10, the whitecap coverage, and the coverage with streaks made of sea foam, bubbles, and spray. Their findings are that CD10 increases with wind up to U10 ≈ 40 m s−1 and then decreases with further wind increase. The drag coefficient at U10 ≈ 60 m s−1 is CD10 = 0.0007 ± 0.0003. The whitecap coverage is low, about 4%, indicating a low level of the dominant-scale breakers. In contrast, the streak coverage increases nearly to 100% at U10 ≥ 40 m s−1. It is suggested that the streaks form a layer of water–air emulsion at the air–sea interface, and such a layer acts as a slip layer for the wind (Powell et al. 2003; Holthuijsen et al. 2012). Our LES results (Figs. 2c,d) also show that those sea states that yield CD10 ≈ 0.0007 at U10 = 59 m s−1 must have a low level of the dominant-scale breakers (i.e., fb ≲ 1) as well as a slip-unresolved surface (i.e., ). According to our results, a decrease of CD10 with a wind increase as seen in Holthuijsen et al. (2012) occurs only if at least one of the sea state parameters decreases.
The drag coefficient CD10 is sensitive to the breaker factor fb in the sense that increasing fb only by 6 times is enough to make a CD10 variation comparable to that in the observational data. The sensitivity of CD10 on fb can be seen from the slopes of the isocontours in Figs. 2b, 2c, and 2d; that is, the isocontours become more horizontal when CD10 is more sensitive to change in fb. The sensitivity significantly increases (i.e., CD10 increases more rapidly with increasing fb) as the wind speed increases (from Fig. 2b to 2c) and as the dominant scale of the breakers λo increases (from Fig. 2c to 2d). This result emphasizes the importance of accurate observations of the breaking wave distributions in high winds.
b. The importance of the dominant-scale breakers in momentum exchange
The ratio of the stress supported by the dominant-scale breakers to the net wind stress is shown by color in Figs. 2b, 2c, and 2d. Regardless of the dominant-scale λo, the dominant-scale breakers support a significant percentage (viz., 30%–100%) of the net drag at U10 = 59 m s−1 at those sea states that yield a CD10 close to the lower bound of the observational CD10 range (Figs. 2c,d). For those conditions that yield higher values of CD10, the importance of the dominant breakers varies. However, when the breaker factor fb is more than 1 and λo is longer (Fig. 2d), the dominant-scale breakers support significant drag.
At U10 = 17 m s−1, the importance of the dominant-scale breakers is reduced (Fig. 2b). When fb is about 1, the dominant-scale breakers are not important. However, their importance quickly increases with increasing fb. Considering the fact that the natural variability of Λ can be an order of magnitude at this wind speed (Kleiss and Melville 2010), there may be plenty of sea states where the dominant-scale breakers are important.
c. Mean wind profile
The normalized mean wind shear is presented to show the shape of the mean wind profile and its dependence on the wind and sea state (Figs. 3a,b). In particular, the results shown are obtained at = (3.5, 0.0006) and (1, 0.0022) for three combinations of the wind and the dominant scale: namely, (U10, λo) = (17 m s−1, 29.22 m), (59 m s−1, 29.22 m), and (59 m s−1, 256.46 m). These conditions of U10 and λo are the same as those in Figs. 2b, 2c, and 2d. The two cases of = (3.5, 0.0006) and (1, 0.0022) are representative of the sea states in which the impact of the dominant-scale breakers is significant and insignificant, respectively.
The von Kármán constant κ obtained from the best fit is always 0.37; this value is within its typical range 0.35 ≤ κ ≤ 0.44 (Andreas 2009). The results of ϕm show insignificant dependence on U10 and λo at a given . In contrast, parameters significantly influence ϕm. For reference, the normalized mean wind profile at U10 = 59 m s−1 and λo = 256.46 m is shown in Figs. 3c and 3d. When the dominant-scale breakers support less than 60% of the net drag (i.e., Figs. 3b,d), the mean wind stays more or less logarithmic even in the roughness sublayer (viz., z/ao ≲ 8). The impact of the dominant-scale breakers on the mean wind profile becomes discernible when the dominant-scale breakers support more than 60% of the net drag (i.e., Figs. 3a,c). Their main impact appears below z/ao ≈ 0.5 and is to reduce the mean wind shear and increase the mean wind speed compared to the corresponding log profile. Such a result is similar to wind profiles over other types of tall roughness (e.g., Britter and Hanna 2003).
There is a notable difference between the mean wind profile over the (idealized) open ocean in this study and that over a very short (viz., λo ≲ 1 m) and young wave field typical in laboratory wind-wave tanks (Suzuki et al. 2013). Over the laboratory water surface, ϕm around z/ao = 1 is significantly larger than 1. This is because the wave breaking of the very short dominant waves in the laboratory is much more frequent than the breaking of the dominant waves in the open ocean. In the laboratory, the regions of separated flow induced by the dominant-scale breakers occupy the majority of the near-surface layer, and the mean wind shear around z/ao = 1 is dominated by the high wind shear that develops along the separated flow regions (Suzuki et al. 2013). In contrast, over the open ocean, much of the near-surface layer is outside the separated flow regions made by the dominant-scale breakers, and the statistics of the mean wind profile is dominated by the flow pattern outside the separated flow regions.
d. The impacts of dominant-scale breakers on turbulence characteristics
According to Suzuki et al. (2013), the main types of the turbulence in the boundary layer over a very short and young breaker field are quasi-streamwise vortices (i.e., regular shear turbulence) and breaker-induced wake turbulence, and the statistical properties of the near-surface turbulence result from a mixture of their turbulence characteristics. The same picture holds true at the (idealized) open-ocean conditions in the current study. However, there are significant differences in the turbulence statistics between the two studies because the dominant-scale breakers at the open ocean occur infrequently compared to the laboratory ones. The dominant-scale breakers at the open ocean are usually not located in other breakers’ wakes and, hence, are fully exposed to high wind. As a result, they produce very intense wake turbulence.
Such intense wake turbulence is characterized with large negative (Fig. 4a), large TKE (Fig. 4b), large upward wind burst (Fig. 4c), large upward TKE flux (Fig. 4d), and large ejection (i.e., where and ) (Fig. 4e). Here, the TKE is defined as , and the TKE flux is defined as
There is also a large signature in the stress (Fig. 4f); a strong downward flux is followed by a strong upward flux downstream.
The dominant-scale breakers not only produce low speed wakes but also increase the surface wind outside the wakes; that is, is positive and very large outside the wakes (Fig. 4a). Correspondingly, the mean wind near the surface is higher than the log profile counterpart (Fig. 3c). This increased surface wind will counteract the effect of the reduced wind speed (sheltering effect) occurring inside the wakes.
The wake turbulence is so intense that, although its occurrence is infrequent, it can dominate the statistical characteristics of the turbulence in the roughness sublayer. Figure 5 shows the distinct impacts of the dominant-scale breakers in the normalized variances, the averaged TKE flux normalized with , and . Again, the impacts of the dominant-scale breakers are significant at = (3.5, 0.0006) and insignificant at = (1, 0.0022). Thus, the effects of the dominant-scale breakers can be seen by comparing the results at = (3.5, 0.0006) to those at = (1, 0.0022). The intense wake turbulence enhances the variances below z/ao ≈ 1, making the near-surface layer about twice as turbulent at = (3.5, 0.0006) than at = (1, 0.0022).
The dominant-scale breakers also make more positive and larger below z/ao ≈ 9 (cf. Figs. 5c and 5d). The quantity indicates the balance between positive and negative . The large positive in Fig. 5c confirms the fact that the dominant-scale breakers cause bursts of upward (i.e., away from the surface; ) airflow where the flow separates from the surface (Fig. 4c). These upward wind bursts may have an important implication on suspension and dispersion of sea spray, aerosols, and other tracers. The generation of such tracers occurs largely in wave breaking events, where intense tearing of sharp wave crests, splashing of water, mixing of air and water, and bursting of bubbles take place. Once generated, they may be selectively transported away from the surface due to the upward bursts of wind.
The TKE flux is another quantity that is heavily influenced by the dominant-scale breakers (cf. Figs. 5e and 5f). For = (1, 0.0022), the TKE flux is negative adjacent to the surface and becomes positive and stays relatively constant away from the surface. This is similar to the TKE flux in other flows (Ikeda and Durbin 2007; Lee and Sung 2007). In contrast, the dominant-scale breakers cause large positive TKE fluxes at the points of flow separation (Fig. 4d), make the overall TKE flux positive even adjacent to the surface, and enhance it nearly 3 times at heights below z/ao = 6 (Fig. 5e). These local TKE fluxes induced by the dominant-scale breakers transport part of the wake TKE away from the surface (Figs. 4b,d and 5e).
The normalized turbulence statistics in Figs. 5a, 5c, and 5e are dominated by the wake turbulence properties and reflect the normalized intensity of the wake turbulence. In all of these figures, the dotted lines (i.e., the normalized turbulence quantities at U10 = 17 m s−1) are always less than the dashed and solid lines (i.e., the normalized turbulence quantities at U10 = 59 m s−1). Therefore, the normalized intensity of the wake turbulence increases with increasing wind.
e. TKE budget
In a statistically steady state, the TKE satisfies (Suzuki et al. 2013)
where and is the work done on the turbulence by the breaker-induced pressure perturbation and represents the production of wake turbulence. It satisfies . The terms on the right-hand side of Eq. (10) are called the TKE transport, shear production, wake production, and viscous dissipation, respectively. The balance of these terms at U10 = 59 m s−1 is shown in Figs. 6a and 6b.
Figure 6a shows four major effects of the dominant-scale breakers in contrast to Fig. 6b. First, the dominant-scale breakers produce intense wake turbulence. The wake production is far more important than the regular shear production near the surface. Second, they reduce the shear production. This is because they reduce both the mean wind shear and the Reynolds stress near the surface. The Reynolds stress inevitably reduces with decreasing height near the surface as part of the momentum being transferred down by the Reynolds stress is absorbed by the breakers via the breaker-induced pressure perturbation (form drag) (Suzuki et al. 2013). Third, they induce upward TKE fluxes described in section 3d. These TKE fluxes take away a large amount of wake TKE at 0 < z/ao < 0.5 and deposit it at z/ao > 0.5. The TKE transport term can be as significant as other terms. Fourth, the TKE dissipation at the resolved heights can be enhanced by the dominant-scale breakers. This is because the increase in the wake production can exceed the reduction in the shear production.
There is a Reynolds-averaged Navier–Stokes (RANS) parameterization (Hara and Belcher 2004; Kukulka et al. 2007; Kukulka and Hara 2008a,b) of the TKE dissipation that assumes that the TKE dissipation is proportional to the turbulent stress raised by the exponent of . In Suzuki et al. (2013), it was pointed out that this RANS parameterization significantly underestimates the TKE dissipation in the roughness sublayer over a very short and young breaker field. Such underestimation of the TKE dissipation causes an overestimation of the drag coefficient in the RANS. Comparison of Figs. 6c and 6d suggests the same conclusion in open-ocean conditions when the breaking wave effect is significant [i.e., = (3.5, 0.0006)].
Here, we will discuss the effect of breaking waves on the drag coefficient in detail. Consider the atmospheric near-surface layer between the water surface and some height h inside the log layer. The total stress (divided by the density) is constant at in this layer. The drag coefficient CDh with respect to the mean wind speed Uh at z = h is defined as CDh = (U*/Uh)2. The ratio Uh/U* is equal to a normalized energy flux (density) , which is the rate of work done at the layer top by the stress on the mean wind normalized by . This normalized energy flux into the mean wind through the layer top is downward; that is, it causes a gain in the mean wind energy of the near-surface layer. According to the definition, the drag coefficient CDh is smaller when this normalized energy input to the near-surface layer is larger.
In a statistically steady state, the energy input to the near-surface layer is larger when the near-surface layer fluxes out or dissipates more energy. This can be shown by the kinetic energy budget integrated over the near-surface layer and normalized with (Suzuki et al. 2013): namely,
where Pwave is the energy transfer to the waves (wave production) via the work done on the waves by the wind forcing. The lhs is the normalized energy input to the layer, and the terms on the rhs are the normalized energy outputs from the layer: namely, the outflux of TKE through the layer top (first term), the TKE dissipation over the layer (second term), and the energy transfer to waves (third term). In summary, Eq. (11) shows that, for a given wind stress, the reference wind is higher or the drag coefficient is lower when the near-surface layer fluxes out or dissipates more energy.
At all conditions considered in this study, the TKE dissipation is the largest term on the rhs of Eq. (11). As shown in section 3e, when the effect of the dominant-scale breakers increases, the wake production becomes large enough to replace the reduction of the shear production. Therefore, the net (i.e., the sum of the shear and wake) TKE production stays large and keeps the TKE dissipation large. As a result, the drag coefficient CDh remains relatively small.
Figure 7 is a diagram depicting the energy flux balance shown in Eq. (12). At the layer top, the mean wind gains energy. Part of that energy converts into large- and small-scale eddies via the regular shear turbulence processes. When the mean wind forces a wave via the pressure perturbation induced around the wave and wake turbulence is produced as a result, the mean wind transfers its energy to the wave and the wake turbulence. A similar energy transfer occurs when a gust (large eddy) forces a wave, but this time the gust TKE converts into the wave energy and the wake TKE. The scale of the wake turbulence is comparable to the wave height or less. Hence, generation of wake turbulence transfers energy from very large scales (mean wind and gusts) to relatively small scales and shortcuts the regular energy cascade.
Although presenting the complete picture of the energy flux partition [Eq. (12)] is beyond the scope of our current study, let us consider how these terms may depend on the wind. A conceptual schematic of the wind dependence is shown in Fig. 8. At very low winds, waves are not forced strongly. Hence, the impact of waves is small, and the energy flux partition is similar to that of a flow over a flat wall; that is, the wave and wake productions are small, and most of the energy influx at the layer top goes into TKE via the shear production (similar pattern to Fig. 6b). As the wind increases, the forcing on the waves increases and so do the wave and wake productions. At the same time, however, the waves would impede the normalized shear production as the momentum transfer to the waves reduces the normalized Reynolds stress (and usually the normalized mean wind shear) near the surface, similarly to the impact of the momentum transfer to breakers seen in section 3e. The increase of the drag coefficient with increasing wind typical in observations at low-to-moderate winds implies that this decrease in the normalized shear production exceeds the increase in the normalized wave and wake productions. This is probably related to the fact that the impact of the waves at low-to-moderate winds is still confined very close to the water surface. As the normalized shear production very close to the surface is very large without the presence of waves [viz., it is 1/(kz)], the diminishing of the shear production there by the waves has a large impact. In contrast, since the wind very close to the surface is small the wake production tends to be small and cannot replace the reduction in the shear production.
Our LES results at U10 = 17 m s−1 and λo = 29.22 m show that the contribution of the dominant-scale breakers to the normalized wake production [second term on the rhs of Eq. (12)] is only 2% of the net energy flux [viz., the lhs of Eq. (12)] at sea state = (1, 0.0022) and 30% at = (3.5, 0.0006). The contribution of the dominant-scale breakers to the wave production [third term on the rhs of Eq. (12)] for the same wind-sea conditions is 2% and 18%, respectively.
At high winds, the normalized wake production significantly increases. In contrast, further decrease in the normalized shear production is not as drastic probably because the shear production near the surface is already impeded fully as seen in Suzuki et al. (2013). In our LES results at U10 = 59 m s−1 and λo = 256.46 m, the contribution of the dominant-scale breakers to the normalized wake production is 10% of the net energy flux at sea state = (1, 0.0022) and 50% at = (3.5, 0.0006). Their contribution to the normalized wave production for the same wind-sea conditions are 4% and 16%, respectively. Without the high wake production in Fig. 8, the drag coefficient would be much larger.
The impacts of the dominant-scale breakers on the air–sea momentum flux and boundary layer turbulence are investigated using LES in open-ocean conditions in moderate-to-high winds (viz., 17 ≤ U10 ≤ 59 m s−1). At many sea states likely to occur in such conditions, the dominant-scale breakers are found to support a significant percentage of the total wind stress. This is an interesting result since it is often suggested that the dominant-scale breakers are unimportant at the open ocean, even in high winds, except for very short fetch conditions (Makin and Kudryavtsev 2002; Kudryavtsev and Makin 2007; Mueller and Veron 2009).
When the sea state parameters are smaller and are fixed, CD10 increases with the wind until it reaches an upper-bound value. Then it saturates with further wind increases. A decrease of CD10 with increasing wind as seen in Powell et al. (2003) and Holthuijsen et al. (2012) occurs only if the unresolved roughness and/or the dominant-scale breakers reduce with increasing wind.
The dominant-scale breakers induce airflow separation. Outside the separated flow regions, the local surface wind is higher and the local wind shear is lower than the logarithmic wind profile having the same CD10 and U10. As the dominant-scale breakers occur infrequently, much of the near-surface layer is outside the separated flow regions. Hence, the flow pattern there determines the mean wind profile. The effect of the increased wind speed outside the separated flow regions will counteract the effect of the reduced wind speed (sheltering effect) occurring inside the separated flow regions.
In general, CD10 is sensitive to the amount of the dominant-scale breakers. For example, increasing their amount by six times results in a CD10 variation that roughly spans the range of CD10 variation typically observed in the open ocean. The sensitivity increases with increasing wind or increasing breaker-scale λo. This result emphasizes the importance of accurate observations of the breaking wave distributions in moderate-to-high winds.
Because they are infrequent and fully exposed to high wind without locating in other breakers’ wakes, the dominant-scale breakers produce very intense wake turbulence. Especially the ejections and upward wind bursts induced by the breakers can transport significant amount of TKE away from the surface. Moreover, as a large release of sea spray, aerosols, and other tracers to the air takes place mostly near the breaking crests, these turbulent motions may have important implications to suspension and upward transport of these tracers.
The dominant-scale breakers impede the regular shear production of TKE and, instead, produce wake turbulence. The normalized wake production increases with wind, and it can dominate the shear production at many sea states in high winds. At the hurricane strength wind of U10 = 59 m s−1, the wake production by the dominant-scale breakers reaches 10%–50% of the total energy flux into the surface layer below 10-m height. Hence, ignoring the wake production may cause a severe overestimation of CD10.
This work was supported by the U.S. National Science Foundation (Grant OCE0824906). We used the computational resources at National Center for Atmospheric Research (NCAR).
Current affiliation: Brown University, Providence, Rhode Island.