TMCnet News

Large-Eddy Simulation of Marine Atmospheric Boundary Layers above a Spectrum of Moving Waves [Journal of the Atmospheric Sciences]
[November 07, 2014]

Large-Eddy Simulation of Marine Atmospheric Boundary Layers above a Spectrum of Moving Waves [Journal of the Atmospheric Sciences]


(Journal of the Atmospheric Sciences Via Acquire Media NewsEdge) ABSTRACT Momentum and scalar transport in the marine atmospheric boundary layer (MABL) is driven by a turbulent mix of winds, buoyancy, and surface gravity waves. To investigate the interaction between these processes, a large-eddy simulation (LES) model is developed with the capability to impose a broadband spectrum of time-varying finite-amplitude surface waves at its lower boundary. The LES model adopts a Boussinesq flow model and integrates the governing equations on a time-varying, surface-fitted, nonorthogonal mesh using cell-centered variables with special attention paid to the solution of the pressure Poisson equation near the wavy boundary. Weakly unstable MABLs are simulated with geostrophic winds increasing from 5 to 25 m s-1 and wave age varying from swell-dominated to wind-wave equilibrium. The simulations illustrate cross-scale coupling as wave-impacted near-surface turbulence transitions into shear-convective rolls with increasing distance from the water. In a regime with swell, low winds, and weak heating, wave-induced vertical velocity and pressure signals are readily observed well above the standard reference height ?a = 10 m. At wind-wave equilibrium, the small-scale wave-induced signals are detectable only near the water surface. Below ?a, a nearly-constant-flux layer is observed where the momentum flux carried by turbulence, form stress, and subgrid-scale motions shifts with varying wave age and distance above the water. The spectral content of the surface form stress is wave-age dependent, especially at low wavenumbers. The LES wind profiles deviate from Monin-Obukhov similarity theory in nonequilibrium wind-wave conditions, and entrainment is greatly enhanced by shear-induced engulfment events.



(ProQuest: ... denotes formulae omitted.) 1. Introduction To go beyond the current statistical description of atmosphere-wave-ocean coupling used in weather and climate models (e.g., Chen et al. 2013; Fan et al. 2012) requires a deeper understanding of the coupling (or interaction) dynamics at time and space scales commensurate with the winds, waves, and currents that exchange momentum and scalar fluxes in the lower atmosphere and upper ocean (Sullivan and McWilliams 2010). Numerous questions remain unsettled about wind-wave coupling processes in the marine atmospheric boundary layer (MABL); for example, the vertical extent of the waveimpacted (or wave induced) layer above the ocean surface, the role of swell, the validity of Monin-Obukhov (MO) similarity theory to predict surface fluxes, the correlation between winds and waves for varying wave state, and how wave-influenced surface wind and pressure fields blend into the bulk of the MABL, to mention just a few. Also there is still vigorous debate as to the mechanisms leading to wave growth, the role of critical layers (Miles 1957; Lighthill 1962; Belcher and Hunt 1998; Hristov et al. 2003; Sajjadi et al. 2014), and the drag of the sea surface at high winds (Bell et al. 2012). Numerous field campaigns have provided important insights (e.g., Edson et al. 2007; Högström et al. 2009; Grare et al. 2013; Högström et al. 2013) but still lack some of the key observations (e.g., pressure-wave slope correlations), as they are forced to cope with the inherent unsteady motion of the wavy interface, the nonstationarity of windwave states, and the wide range of time and space scales involved [e.g., millimeters (spray) to kilometers (remotely generated swell)].

Turbulence-resolving direct numerical simulations (e.g., Sullivan et al. 2000; Yang and Shen 2010; Druzhinin et al. 2012; Richter and Sullivan 2013) and large-eddy simulations (Suzuki et al. 2013; Yang et al. 2013) of idealized wavy Couette flows along with recent laboratory measurements (e.g., Donelan et al. 2004; Veron et al. 2007) have added to the understanding of wind-wave couplings. For example, this body of work highlights the importance of wave age in generating wave-induced motions and the impact of wave breaking, flow separation, and wind speed on surface drag. However, largescale turbulence generated by shear and convection, three-dimensional boundary layer winds generated by rotation, and a broadband spectrum of surface waves, which are fundamental to the MABL, are of course missing in these small-scale studies.


High Reynolds number large-eddy simulation (LES) has the potential to shed light on wind-wave coupling processes in the MABL, but its application is relatively unexplored because of the complexity of including a time-evolving wave field with multiple modes and the high computational expense required to span a wide range of scales. The topic of the present paper is the latter (i.e., to develop and apply LES to a realistic stratified MABLwhere the lower boundary is a resolved spectrum of time-dependent surface gravity waves). The computational method described here allows for nearly arbitrary 3D wave fields [i.e., the sea surface elevation h 5 h(x, y, t) is a surface boundary condition; h(x, y, t) is assumed to be a single-valued function]. The spatial scales of the resolved turbulence and waves considered are on the order of meters up to the height of the MABL zi ; O(500 m) or more. Covering this broad range of scales requires a large number of mesh points and significant computational power. In the current problem posing, the surface waves are externally imposed based on existing empirical wave spectra with the assumption of a linear-mode random-phase model for the wave amplitudes. The missing phase information can in principle be obtained by using advanced observational techniques of the sea surface (Romero and Melville 2010). The present work continues past efforts that coupled the MABL with an idealized single monochromatic wave (Sullivan et al. 2008; Nilsson et al. 2012), but it is a significant advance over these earlier developments.

The outline of the paper is as follows: section 2 is an introduction to the LES equations appropriate for a high Reynolds number atmospheric boundary layer and their transformation to time-dependent curvilinear coordinates, section 3 outlines the numerical method, section 4 describes results for a weakly unstable MABL with varying geostrophic winds, and section 5 provides a summary of the current findings.

2. Governing equations In the description of the model equations, given below, the following notation is used: ui5(u, y, w) denotes the Cartesian velocity components, u is virtual potential temperature, p* is the pressure variable normalized by density r, and e is the subgrid-scale (SGS) energy. Flow variables (ui, u, p*, e) are interpreted as LES spatially filtered quantities.

a. Atmospheric LES model equations The set of spatially filtered LES equations applicable to turbulent flow in the atmospheric boundary layer under the Boussinesq assumption is (e.g., Moeng 1984; McWilliams et al. 1999) ... (1a) ... (1b) ... (1c) ... (1d) In the above, Vi 5 (0, 0, VsinF), with Earth's rotation and latitude denoted V and F, respectively; the buoyancy parameter b 5 g/uo, where the gravitational acceleration is g and the reference (still air) virtual potential temperature is uo; and the large-scale external pressure gradients normalized by density ?P/?xi drive the boundary layer winds. A Boussinesq flow model requires that p*, the pressure variable, be determined from an elliptic Poisson pressure equation to enforce mass conservation (see section 3d). Ultimately, the SGS momentum and temperature (or scalar) fluxes (tij, tiu), respectively, require modeling in the interior of the flow and at the lower boundary (see section 3f) to close the system of equations. In the SGS turbulence kinetic energy (TKE) equation [see (1d)], the time evolution of e depends on right-hand-side terms (in order): advection by the resolved field, energy transfer between resolved and SGS motions where the resolved scale strain rate Sij 5 (1/2)(?ui/?xj 1 ?uj/?xi), SGS buoyancy flux, Laplacian diffusion with turbulent eddy viscosity nt, and viscous dissipation E.Molecular diffusion terms are neglected assuming that the molecular Reynolds number is high. A discussion of the modeling of the SGS fluxes is postponed until section 2d.

b. Coordinate transformation We adapt our LES model with a flat boundary (Sullivan and Patton 2011) to the situation with a threedimensional time-dependent lower boundary with shape h5h(x, y, t) by applying a transformation to the physical space coordinates xi [ (x, y, z) that maps them onto computational coordinates ji [ (j, h, z). The computational mesh in physical space is surface following, nonorthogonal, and time varying. Vertical gridlines are held fixed at a particular (x, y) location on the surface but translate vertically as a function of time t. A transformation that obeys these constraints and maps the physical domain to a flat computational domain xi0ji is the rule: ... (2) The differential metrics ?xi/?jj and ?ji/?xj, which are needed in formulating the LES model in curvilinear coordinates, are connected through the mapping transformation. This can be done for general transformations (e.g., Anderson et al. 1984, 252-253), but here we take advantage of the simplified terrain-following grid in (2). This leads to the reduced set of nonzero metric relationships: ... (3) where J is the Jacobian and subscripts denote partial differentiation. The time dependence of the mapping appears in (3), where ?z/?t 5 zt is the grid speed (i.e., the vertical velocity of individual grid points). We define the rule in (2) that maps physical space onto flat computational space at any time by the following prescription: ... (4) where ZL is the top of the computational domain and - controls how rapidly the horizontal gridlines become level surfaces in physical space. Typically we use - 5 3, which vertically squeezes (stretches) near-surface grid cells at the wave crests (troughs) by a modest amount, approximately 20%. For - 5 1, (4) reduces to z 5 (z 2 h)/(1 2 h/ZL), which is commonly used to map stationary terrain onto flat computational coordinates (e.g., Henn and Sykes 1999). Equation (4) is a singlevalued mapping that produces smooth continuous metric derivatives (?ji/?xj, ?ji/?t) in the interior of the computational domain depending on the smoothness of the boundary shape derivatives (ht, hx, hy).

c. LES equations in time-dependent curvilinear coordinates The transformation of the LES equations to curvilinear coordinates is obtained by applying the chain rule for differentiation and making frequent use of the metric identity (e.g., Anderson et al. 1984, p. 254): ... (5) To cast the equations in fully conservative form, however, requires care in treating the time derivative. For example, the time derivative of density transforms as ... (6) provided the grid movement obeys ... (7) Equation (7) is a restatement of the so-called geometric (or space) conservation law (GCL) and should be satisfied by the numerical discretization in order for the scheme to be conservative (Thomas and Lombard 1979; Demirdzic and Peric 1990).1 For our surface following grid the GCL simply reduces to (8b) since ?ji/?t 5 (0, 0, 2ztJ). Section 3c describes how (7) is used to track the motion of the grid.

The complete set of LES equations in computational coordinates under the time-dependent surface-following transformation (2) and (3) is then ... (8a) ... (8b) ... (8c) ... (8d) ... (8e) ... (8f) where (8a) is the mass conservation (continuity) equation, (8c) is the momentum transport equation, (8d) is the transport equation for potential temperature, (8e) is the subgrid-scale energy transport equation, and (8f) is the pressure Poisson equation. The right-hand sides of (8c)-(8e) model physical processes in the atmospheric boundary layer, namely, pressure gradients, Coriolis rotation, divergence of subgrid-scale fluxes, buoyancy, and, in the case of the SGS e equation, also turbulent diffusion and viscous dissipation. For completeness these terms are gathered here: ... (9a) ... (9b) ... (9c) Momentum and scalar advection are compactly written in strong flux-conservation form using the volume flux or ''contravariant flux'' velocity: ... (10) whereUi5(U, V,W) are normal to a surface of constant ji (see Fig. 1).

The time dependence of the grid modifies the LES equations: the Jacobian appears inside the time tendency of each transport equation and, as anticipated, the total vertical flux of a variable c depends on the difference between the contravariant flux velocity and the grid speed (W 2 zt)c. Examination of (8) shows that if the velocity and scalar fields are set to constant values then the left-hand sides of (8c)-(8e) reduce to (8b). Hence, the numerical method needs to satisfy the reduced form of the GCL discretely in order to prevent artificial sources and sinks from developing in the computational domain (see section 3c).

d. Parameterization of SGS motions Formally, the LES equations are obtained by applying a spatial filter to the full governing equations (e.g., Wyngaard 2004). The filter is assumed to be homogeneous in space and time, and, as a result, the operation order of differentiation and filtering can be interchanged. For boundary layer flows, this assumption is violated and a commutation error occurs near walls (Ghosal and Moin 1995), but the error is almost always ignored in LES implementations as the inaccuracies in SGS wall modeling are viewed as of greater importance (Sullivan et al. 2003). The parameterization problemis evenmore pronounced with a multimode moving lower boundary because of filtered products of fluctuating metric coefficients and velocity (Chalikov 1978), and extra stresslike terms arise from filtering the boundary shape (Nakayama et al. 2009). In the present application, the flow is assumed to be laterally periodic and our wavefollowing grid transformation results in x 5 j and y 5 h. Thus, we adopt the simplest approach and employ the filtering operation and standard SGS models as in our flat LES code but account for the varying filter width. The LES equations are closed using the SGS parameterizations outlined by Deardorff(1972) and analyzed by Moeng and Wyngaard (1988). Also, the solutions are explicitly filtered (i.e., dealiased) in j-h planes at the end of each time step. This parameterization utilizes an eddy viscosity for SGS momentum and temperature fluxes, and the Lilly-Kolmogorov model for viscous dissipation: ... (11) The turbulent eddy viscosity nt 5CkD ffiffiffipe and turbulent diffusivity nh 5 (1 1 2'/D)nt evolve under the action of the transport equation for SGS TKE [see (8e)]. The constants (Ck, C«) 5 (0.1, 0.93) and corrections to the stability corrected length scale ' are found in Moeng (1984) and Moeng and Wyngaard (1988). The filter length scale D is computed from the averaging volume D3 5 (3/2)2(DjDhDz/J) with the 3/ 2 factor accounting for dealiasing. Because of movement and the stretching- squeezing of vertical gridlines, D varies with position and time. To reduce the reliance on the SGS model, we use fine resolution near the surface (see section 4a). Since the large wave motions are fully resolved, the total TKE at the surface is often large and dominated by resolved variances, further reducing the impact of the SGS model compared to a situation with flat stationary walls.

e. Wave field prescription To complete our LES model for the MABL, we need to prescribe the height of the surface wave field h(x, y, t) in physical space over a typical horizontal domain of the LES. In the absence of a full description of the wave field kinematics and phase relationships, we adopt simplifications that allow us to use existing information. Time and space maps of h(x, y, t) are constructed based on typical empirical fits of measured two-dimensional wave spectra: ... (12) where the amplitude S(k) and directional D(k, f) spectra depend on the magnitude of the horizontal wavenumber k5jkj5jkxi 1ky jj and wave direction f. We choose the functional forms for S(k) and D(k, f) proposed by Donelan et al. (1985) [see also p. 187 of Komen et al. (1994)], which depend on bulk environmental parameters, namely, the reference average surface wind speed Ua at za 5 10 m, the phase speed of the peak in the wave height spectrum Cp, their ratio (the wave age) Cp/Ua, and the mean wave propagation direction hfi. To use these measured spectra, which are expressed in terms of radial frequency v, in the LES they are transformed into wavenumber space using the linear dispersion relation jkj 5 v2/g and properly weighted such that the resulting wavenumber spectra are variance preserving (Phillips 1977, p. 105). In physical space, the rule for constructing the synthetic wave field h(x, y, t) is as a sum of linear plane waves, each with random phase u: ... (13) where the amplitude A2(k) 5 S(k)D(k, f)dkxdky. The random phase model in (13) is computationally advantageous since only an initial wave state h(x, y, to) is required and it also easily adapts to the horizontal dimensions of the LES domain. Two-dimensional fast Fourier transforms (FFTs) are used to compute future values of h(x, y, t) as well as the time and space derivatives of the boundary shape (ht, hx, hy), which are needed to construct the moving mesh and the surface boundary conditions.

f. Velocity boundary condition Both the wind and wave fields are assumed to be spatially periodic in computational j-h planes and at z 5 0 a boundary condition on the winds is imposed consistent with the surface wave field. The total time rate of change of the wave height is ... (14) where (uo, yo, wo) are the water motions. Initially, we assume these water motions are dominated by the wave orbital velocities [i.e., there is no net current; see discussion by Banner andMelville (1976) and Banner (1990)]. The contravariant velocity component normal to the water surface is given by (10): ... (15) Matching (14) and (15) and requiring no flow normal to (across) the wavy boundary leads to ... (16) Equation (16) effectively implies that spray and bubbles do not pass through the interface in our problem. At the upper boundary of the computational domain, the horizontal gridlines are level surfaces in physical space with zj 5 zh 5 0. Then ... (17) so that the resolved vertical flux of momentum and scalars is zero across the top boundary. The specification of the surface fluxes appropriate for a high Reynolds number rough-wall LES model taking into account the wave motion is described in section 3f.

3. Numerical method a. Spatial discretization and variable layout The LES set of (8a)-(8e) is discretized using the welldeveloped techniques in our flat bottom boundary codes (Moeng 1984; Sullivan et al. 1996). Spatial derivatives in the (j, h) computational coordinates are estimated using pseudospectral approximations while centered secondorder- accurate finite-difference approximations are employed for z derivatives. To be compliant with (8b), the spatial differencing evaluates the advective terms in (8c)-(8e) in flux form as opposed to the rotational form used by Moeng (1984) or skew-symmetric form used by Sullivan et al. (2000). The variable layout in the computational mesh uses a collocated arrangement for the fundamental solution variables (ui, p*, u, e), which is advantageous as it results in a compact differencing stencil and is also applicable to meshes with large bends in the coordinate lines (Zang et al. 1994; Sullivan et al. 2000). The fundamental difficulty of tightly coupling the velocity and pressure in a collocated mesh is circumvented by locating the contravariant flux velocities at cell faces mimicking the usual arrangement in a staggered mesh as shown in Fig. 1.

b. Time advancement Time integration is a fully explicit third-order Runge- Kutta (RK3) scheme (Spalart et al. 1991; Sullivan et al. 1996, 2008) that uses dynamic time stepping with a fixed Courant-Fredrichs-Lewy (CFL) number and employs a fractional step method to enforce the divergence-free condition. The general rule for advancing a cellcentered Cartesian velocity variable to a new time level n over a time step Dt is ... (18) where the pressure gradient is written in conservative form. The intermediate velocity ... (19) is the discrete sum of the full right-hand-side Qi/J from the previous time (or stage) step (n 2 2) and the partial right-hand-side qi/J from the current time (n 2 1) minus the pressure contributions; an and bn are weights associated with the RK3 scheme (Sullivan et al. 2000). The time integration rule for scalar fields is identical to (18) and (19) in the absence of pressure gradients. The divergence-free condition is satisfied at every RK3 stage, as opposed to the method advocated by Le and Moin (1991) (see also Zheng and Petzold 2006), which projects the velocity onto a divergence-free field only at the final stage. The latter scheme is computationally more efficient since it reduces pressure iterations over a full time step, but both flux and skew-symmetric forms of the advection term in our scalar equations are based on satisfying (8a) at all stages of the RK3 time advancement. Attempts to use the method of Le and Moin (1991) resulted in numerical instabilities.

c. Geometric conservation law The geometric conservation law [(8b)] is satisfied using the same time-advancement scheme as for velocity and scalar fields [see (18)]. Given an externally imposed wave field, h(x, y, t) is known at current and future time steps (n 2 1, n) and we can then construct the distribution of vertical grid points at ''W points,'' [viz., z(x, y)(n21,n) according to the map of (4)]. Jacobians J(n21,n) are then naturally built from the transformation rule in (3). Since we choose to specify the grid point locations at future time steps, we then find matching grid speeds zt(x, y) so that (8b) is obeyed discretely. The companion to (18), appropriately inverted, for determining grid speeds is the semidiscrete relationship: ... (20) This first-order equation is integrated vertically from the surface to the top of the domain with the boundary condition zt 5 ht at z 5 0. We emphasize that the upper boundary of our computational domain is far from the lower surface and the wave-following gridlines in the interior of the domain quickly asymptotically approach flat level surfaces with increasing z. As a result, zt at a particular x-y grid point must vary with distance from the wave surface. This is different than the scheme proposed by Chalikov (1998) in which the grid speeds at all z are equal to the vertical motion of the wave field (i.e., zt is constant with height). The latter scheme also implies that the computational mesh at any z, including the upper boundary, mimics the shape of the lower surface, which complicates posing upper boundary conditions. There is a subtlety in the implementation of (20). Since our time stepping is explicit, the value of zt diagnosed from (20) is at time level n 2 1 but uses gridpoint locations at time level n. In other words, if the wave field is imposed, diagnosing the grid speed at time n requires knowledge about the wave shape at a future time n 1 1. The end result is that the grid must be stored at three time levels in the LES code. Information about the wave field at each RK3 stage is used in two ways in the grid generation: h(x, y, t) is used to position the grid nodes in the mesh and ht(x, y, t) is used to determine how rapidly the grid nodes in the mesh move to their new locations.

d. Pressure Poisson equation The heart of the numerical scheme is the formulation and solution of a pressure Poisson equation that results in a flow field consistent with our incompressible Boussinesq flow model. As is standard practice, we develop a pressure equation by combining the continuity equation and our time-stepping scheme. Substitution of (18) into (10) leads to the time advance rule for the contravariant velocity: ... (21) where now the pressure gradient is expanded in nonconservative form. Momentum interpolation (Sullivan et al. 2000) of the Cartesian velocity components is used to build ... (22) where (*)jI denotes an interpolated value. As mentioned previously, because of the spatial differencing scheme, only W jn21, located at the upper and lower cell faces, requires interpolation (see section 3a). The divergence of (21) leads to the pressure Poisson equation: ... (23) Our strategy for finding the pressure utilizes an iterative method that avoids directly forming the complicated left-hand side of (23). Instead, examination of (18), (21), and (23) suggests the equivalent stationary iteration scheme to find p*,n: ... (24) where the source term is the divergence-free condition, superscript i denotes the iteration index, and the linear preconditioning operator is ... (25) At each iteration, we update Uk using the latest p*,i field according to (21) and then compute the new right-hand side of (24). The operator L is an easily invertible diagonal approximation of the left-hand side of (23) with the spatially averaged Jacobian ... (26) Equation (24) is solved using standard methods, namely, 2D Fourier decomposition in j-h planes followed by tridiagonal matrix inversion in the z direction. Inversion of (24) is expensive because its computational steps require global communication in a parallel algorithm (Sullivan and Patton 2011). The iterative solution for pressure in the presence of wavy boundaries increases the overall computational cost of simulations by approximately 50% compared to simulations with flat lower boundaries where the pressure is solved for directly.

e. Pressure boundary conditions at a wavy boundary Compatible velocity and pressure boundary conditions are naturally built into the iteration scheme given by (24). First, the proper surface boundary condition on the diagonal preconditioner L is simply ?p*/?z 5 0, which is the condition applied in a flat-wall LES with a Cartesian grid. The proper pressure boundary condition used in the source term ?Uk(p*)/?jk of (24), however, needs to account for the nonorthogonal mesh [also see Zang et al. (1994), p. 22]. This boundary condition is exposed by examining the update for (U, V,W) in advancing time from n 2 1 to n given by (21): ... (27a) ... (27b) ... (27c) In the interior of the computational domain, the pressure gradient ?p*/?z, which appears in all three equations, is evaluated using standard second-order centered finite-difference formulas, while the pressure gradients ?p*/?j and ?p*/?h are evaluated using pseudospectral differencing. Inspection of (27a) and (27b) shows that updating U and V using a centered formula for ?p*/?z at the first level z 5Dz/2 above the wavy boundary requires the pressure field at z 5 3Dz/2 and at the ghost point z 5 2Dz/2 below the wavy boundary. Along z 5 2Dz/2, U and V are globally connected to an entire plane of ghost-point pressures p*(x, y, 2Dz/2) through the gradients ?p*/?j and ?p*/?h.

To find the proper ghost-point pressures, we use the W equation, its surface boundary condition, and our iteration scheme. At z 5 0, (27c) is rearranged to yield ... (28) where we impose the boundary condition (16). The superscript i in (28) is again the iteration index, and for clarity we drop the superscripts n and n21 denoting the time step. At each iteration with a known field p*(j, h, Dz/2)i, (28) is rearranged to solve for the ghost-point pressure p*(j, h, 2Dz/2)i with the right-hand side (centered on z 5 0) lagged from the previous iteration. Subsequently, (27a) and (27b) are updated at all interior nodes and a pressure iteration is completed by sweeping through (24) with a new right-hand side. The term W (z50), which appears in (28), is constant during pressure iterations and is obtained from interior-node interpolations.

f. Surface fluxes In high winds, the ocean surface is assumed to be in a fully rough regime with negligible contributions from the thin molecular viscous sublayer (Donelan 1998). Here, we follow the approach of second-order closure modeling for turbulent flow over hills and waves and assume that the water surface is fully rough so that the surface fluxes of momentum and scalars can be represented in terms of law-of-the-wall relationships (e.g., Gent and Taylor 1976; Makin et al. 1995; Belcher and Hunt 1998). This approach is less certain in LES modeling where surface undulations (or waves) are partially resolved; that is, the velocity and pressure fields induced by the large energetic components of the wave field are resolved and the smaller-scale, less-energetic waves are unresolved (e.g., Nakayama et al. 2004; Yang et al. 2013). We simply parameterize the effects of the unresolved waves (i.e., those below the cutoffscale of the LES), using a bulk zo roughness. Then the bulk aerodynamic formulas are applied point by point along the wave boundary with this imposed zo roughness to compute the surface momentum and scalar fluxes as described in the appendix.

4. Turbulent winds over waves a. Simulation details The parameter space spanned by the many combinations of large-scale winds, surface heating, and surface waves forcing an MABL is extremely broad and can vary markedly with time as found in the recent High- Resolution Air-Sea Interaction (HiRes) field campaign; see Grare et al. (2013). Over the short 2-week observational period, the HiRes MABL regime featured moderate to high surface-layer winds varying from 5 to 15ms21, with weak stability conditions just left(unstable) and right (stable) of neutral. The high surfacelayer winds generated large waves with significant wave height (Hs . 4.5 m) over several days. And because of the variability in the wind conditions, wave age varied from growing waves (Cp/Ua , 1.2) to swell-dominated conditions (Cp/Ua . 4). This wide and complex observational space can only be sparsely explored with finite computational resources. Thus, for our LES study, we focus on a series of simulations of a weakly unstable MABL, mimicking the forcing and wave conditions in HiRes, with varying geostrophic wind Ug 5 (5, 10, 20, 25)ms21 and constant surface heatingQ*50.01Kms21. The addition of a slight amount of surface heating serves two purposes: it helps to initiate the turbulence and thereby skips the long spinup time of a purely neutral simulation, and second, it allows an investigation of the large-scale MABL motions generated by shear and convection in the presence of surface waves. For the wave spectrum, we choose wind-wave equilibrium conditions Cp/Ua 5 1.2 and Ua 5 15ms21 with the resulting significant wave height Hs56.4 m. This suite of simulations, although idealized, allows us to carry out a process study with wave age varying from near wind-wave equilibrium to swell dominated in the presence of large surface waves. Here wave-age variations are generated solely by varying Ug. For our comparisons, we also generate a solution for stationary waves (i.e., fixed bumps), with Ug 525ms21.

A typical snapshot of the wave field generated by (13) and used in the LES is shown in Fig. 2. Broadly, the image reveals a preference toward long crested waves; that is, waves with aspect ratio long in the direction perpendicular to the winds; a consequence of the directional wave spectrum. A closer examination of the wave field shows rich multiscale variations in both x and y directions. Horizontal 1D spectra of the wave height variance in the x and y directions, given in Fig. 3, confirm the horizontal anisotropy of the wave field at low wavenumber. The spectra exhibit a power law of k25/2 at small scales (high wavenumbers).

To capture the widest possible range of scales in a full MABL we utilize a computational box of size (Lx, Ly, Lz) 5 (3000, 3000, 800)m and discretize the domain using (Nx, Ny, Nz) 5 (1024, 1024, 512) grid points (i.e., 0.5 3 109 grid points). The horizontal grid spacing Dj 5 Dh 5 2.93m and the first vertical level is nominally located 1m above the water surface. A nonuniform vertical grid distribution is generated by smooth algebraic stretching with the spacing ratio between any two adjacent cells Dzi11/Dzi'1.009.Atypical vertical x-z slice of the time-varying surface-fitted grid mesh near the water surface is given in Fig. 2. The image clearly shows how the mesh translates vertically in response to the wide range of wave modes at the lower boundary. We find that over most of the domain the mesh skewness is low as the wave slopes (hx, hy) are typically smaller than 0.35. The initial temperature sounding is constant, u 5 290K up to an inversion height zi 5 300 m; beyond this initial inversion height, u increases linearly at 3 3 1023Km21, and the Coriolis parameter f 5 1024 s21. The unresolved surface waves are simply modeled by a fixed surface roughness zo 5 0.0002 m, noting that the dynamic modeling approach can be used to reduce the sensitivity to the grid resolution (Yang et al. 2013). Random perturbations in velocity and potential temperature are introduced to trigger turbulence. Statistics are obtained by a combination of space and time averages. Because of horizontal periodicity, spatial averages are computed by area averaging over j-h surfaces (i.e., along a surface of constant z). Averaging in a wavefollowing coordinate system allows us to compute statistics down to the water surface (i.e., between and below the wave crests). The resulting vertical profiles are further averaged in time; these averages are indicated by angle brackets. We emphasize that there is no net flow across the wavy lower boundary [i.e.,W(x, y, z 5 0) 5 0], and, hence, no net accumulation of mass in the computational domain similar to flat wall boundary layers. Our time averaging is taken over a time period when the boundary layer is quasi stationary and grows slowly only because of entrainment at the boundary layer inversion. The top of the MABL zi is marked by a steep stable gradient in potential temperature and is estimated using the maximum-gradient method described in Sullivan et al. (1998). Table 1 lists bulk properties of the simulations.

Because of the very high computational expense, these LESs are first started with a flat lower boundary layer (no waves) and run until the turbulence is in near statistical equilibrium. Then the propagating surface waves are ''linearly grown'' to their full height over a period of about 400 s. This strategy significantly reduces the computational cost and also minimizes nonphysical transients in the solutions that can arise by suddenly introducing the surface waves. The simulations with waves are run for more than 100 000 time steps using restart volumes with fully developed turbulence. The iteration count in the pressure Poisson solver is set to 30. Each run consumes more than 106 h running on 2048 computational cores.

b. Flow visualization To gain an overall impression of the flow patterns in the MABL, a small subset of images generated from our simulations is provided in Figs. 4-7. We find vertical velocity w is an excellent diagnostic for visualizing the wide range of scales and wave-MABL couplings captured in our LES. For example, Figs. 4 and 5 illustrate how the vertical velocity field smoothly transitions from a small-scale wave-impacted regime near the water surface into a large-scale streaky structure in the bulk of the MABL. Inspection of the images shows that, at z 5 2.5 m, w tends to closely mimic the long-crested north- south orientation of the streamwise-propagating wave field seen in Fig. 2, this holds for all wave ages. At z 5 9.6 m, the wave-induced pattern in w is visible but gradually blurs with decreasing wave age (i.e., as the background turbulence grows in intensity and the winds and waves approach equilibrium).

At high levels in the MABL, z . 90 m, the dominant flow pattern in w in each simulation switches from the near-surface north-south orientation into a streaky structure mainly aligned in the streamwise (west-east) direction. Animations clearly show the long-lived temporal and large-scale spatial coherence of these structures. Similar shear-convective-generated streaks (or rolls) are present in weakly stratified boundary layers over land in a more unstable range 2zi/L ; 1.5 (e.g., Moeng and Sullivan 1994; Fedorovich et al. 2004) but are pervasive in our simulations of the MABL even with the generally higher winds and weaker surface heating (i.e., 2zi/L , 0.5). Further evidence for shearconvective rolls in the MABL is provided in the crosscutting view (a y-z plane) of Fig. 6 for the simulation very near neutral conditions, Cp/Ua 5 1.1 and 2zi/L 5 0.27. The roll organization over the entire MABL is even stronger for simulations with 2zi/L 5 0.42 and 2.28; this can also be inferred from the w visualizations in Fig. 4. The descending branch of a roll brings high-speed fluid down to the water surface where the descending flow then turns spanwise and converges with a neighboring roll in a layer just above the wavy surface. Finally, eruptions above the wave field liftlow-speed fluid vertically completing a circuit. At za, the ratio of the maximum and minimum values of the horizontal wind Ua(x, y) differs by more than a factor of 2 because of the large-scale circulations. Thus, the lateral heterogeneity induced by the MABL rolls leads to spatial variability in wave age that is larger and more persistent than that due to turbulent wind gusts as shown in Fig. 7. This variability is readily apparent in the w patterns underneath a descending branch of a roll circulation where u(x, y) . hui. The probability density function of wave age shown in Fig. 7 widens considerably in the shear-convective regime (i.e., as the bulk atmospheric stability varies from 2zi/L 5 0.27 to 16.9). Figure 6 shows the large impact of wind shear on MABL entrainment at 2zi/L 5 0.27. Large engulfment events of tropospheric air tilted into the streamwise direction extend deep into the middle MABL (see section 4d and Fig. 13). These engulfment events are frequently observed to be accompanied by intermittent detrainment of MABL air.

c. Winds One of the outstanding questions in MABL dynamics is the interaction between surface waves and the overlying winds, and the vertical extent of the region impacted by surface waves, the so-called wave boundary layer (Sullivan et al. 2008; Sullivan and McWilliams 2010). A separate but related topic is the adequacy of rough-wall Monin-Obukhov (MO) similarity theory, typically applied to homogeneous land surfaces, to predict surface fluxes based on winds measured at a standard reference height. To investigate these questions, we follow a large body of observational studies (e.g., Fairall et al. 2003) and use the LES atmospheric conditions at the reference height za 5 10m to define a nondimensionalizing friction velocity u*,a. First, we compute vertical momentum fluxes uf,a based on time and space averaging in (8c) and (9a) at za: ... (29a) ... (29b) These momentum fluxes include contributions, from leftto right, by resolved turbulence tf, pressure pf, and subgrid-scale turbulence sf terms; (29) includes small but important pressure-wave slope correlations pf to account for the vertical variation of the horizontal gridlines j and h in physical space. The average fluxes in (29) are next aligned with the average winds huai and thus ... (30) serves as our LES estimate of the friction velocity (or wind stress u2 *,a) at za. Note that, in arriving at (30), there is no assumption of a constant-flux layer.

Figures 8 and 9 compare vertical profiles of the nondimensional mean wind speed jhuij/u*,a over the bulk of the boundary layer and also near the water surface for varying wave age. In Fig. 9 we also include a wind profile predicted from (MO) similarity theory (Wyngaard 2010): ... (31) where the stability function c(z/L), given by Large et al. (1994), accounts for the slight surface heating. The vertical coordinate is distance from the water surface z, not from the mean water level. To compare MO theory with LES, we use u*,a defined by (30) and treat zo in (31) as a matching parameter [i.e., we pick zo so that, at the reference height, UMO(za) 5 jhua(za)ij].

Inspection of the bulk wind profiles in Fig. 8 reveals a clear dependence on wave age. In the simulation with wave age Cp/Ua 5 4.2, the winds in the bulk of the boundary layer are faster than for a similar flow above a flat (stationary) boundary. In this nonequilibrium wind-wave regime, the surface-layer winds are accelerated (pushed forward) by fast-moving waves, thereby creating a wave-driven wind. For the wave spectrum in Fig. 3, a large fraction of the wave modes are traveling in the same direction as, but faster than, the local surface-layer winds, roughly all wave modes with l . 16m propagate with speed c . 5ms21, or c/u*,a . 34. The mechanics of the wave-driven wind process are analyzed in a much idealized LES by Sullivan et al. (2008) and Nilsson et al. (2012) where swell is simply modeled as a single monochromatic surface wave. Broadly, we find dynamical similarities in the current LES with a much more realistic full wave spectrum. The key dynamics leading to wave-driven winds is the pressure form stress that acts as a thrust (as opposed to a drag) on the overlying winds at large wave age (see Fig. 10). The overall well-mixed shape of the wind profiles observed in Fig. 8a is a consequence of the relatively light winds forced by a small amount of surface heating; the stability measure 2zi/L 5 16.9 indicates that convection is significant. In our simulations, wave age is varied from 4.2 to 1.1 by increasing Ug from 5 to 25ms21. On the approach to wind-wave equilibrium, the bulk boundary layer winds slow down compared to their counterparts over a flat lower boundary; a consequence of a reversal in the pressure form stress switching sign from positive to negative. With increasing geostrophic winds and surface heating held constant, the wind profiles gradually assume a shape typical of shear-dominated boundary layers (e.g., Moeng and Sullivan 1994) but with higher drag compared to a flat stationary wall, u*,a and zo both increase withUg as shown in Table 1.

A comparison between the LES wind profiles and the MO fit in (31) in the near-surface region further illustrates wave effects (see Fig. 9). As expected, the largest discrepancy between LES and MO wind profiles occurs in the situation with nonequilibrium winds and waves (i.e., light winds above fast-moving waves). The MO fit predicts noticeably faster surface-layer winds below za compared to LES. This discrepancy cannot be simply reconciled by adjustment of the matching height [e.g., choosing za55moverpredicts (underpredicts) the wind speed below (above) za]. Observations reported by Högström et al. (2013) also suggest that the wind profile near the water surface deviates markedly from MO predictions in swell-dominated regimes. The deficiency of MO similarity theory is a consequence of missing dynamics (i.e., upward vertical momentum flux from large surface waves to the atmosphere). Hanley and Belcher (2008) propose a 1D column model that captures some of the effects of wave-driven winds. The agreement between the LES winds and MO fit are remarkably close below za and also in a limited zone above za as the winds and waves approach equilibrium. In this regime, the peak of the spectrum is traveling about 20% faster than the surface winds and then, very broadly, the remainder of the smaller-scale wave modes acts as a slow-moving momentum sink (or drag). The wind profiles displayed in Figs. 8d and 9d further illustrate that even at wind-wave equilibrium conditions flow over moving waves differs from its counterpart over an identical but stationary wave spectrum. In the situation of stationary bumps, the wavy surface is rougher with larger values of u*,a and zo and the vertical height where the wind profile blends into a log-linear profile is higher, near z 5 30 m. In flow over stationary bumps, the winds effectively feel all wave modes as roughness elements. The robustness of the MO fit in (31) clearly needs further evaluation over a broader range of forcings, including more common nonequilibrium wind-wave states (Hanley et al. 2010) using representative wave spectra as discussed here. For example, LES with weak winds opposing waves are markedly different than simulations with wave-driven winds (results not shown).

d. Momentum and scalar fluxes and velocity variances Implicit in the application of MO theory to flow over ocean waves is the assumption that za is located above any wave-induced layer but remains sufficiently close to the water surface so that the change in the average vertical momentum fluxes with height is negligible (i.e., there is a constant-flux layer). In the convective atmospheric boundary layer over land, above the constant-flux layer and below the entrainment zone, the momentum and scalar fluxes decrease almost linearly from their surface values [see p. 215 in Wyngaard (2010)]. The constant-flux-layer assumption is the cornerstone that allows field studies to infer the stress at the water surface using observations at za, and it is also a widely adopted assumption by second-order closure models (e.g., Gent and Taylor 1976; Makin et al. 1995) and in turbulence simulations (e.g., Sullivan et al. 2000; Yang et al. 2013; Suzuki et al. 2013). Here we use LES results to test the validity of the constant-flux layer assumption for the wavy MABL.

Contributions to the u and y momentum equations from the fluxes tf, pf, and sf, essentially the time- and space-averaged right-hand-side terms of (29a) and (29b), are summed and compared in Figs. 10 and 11 for varying wave age over the full depth of the MABL. Remarkably, we find the total flux, computed in wavefollowing coordinates, for the dominant u component is indeed nearly constant over the interval 0 , z , za in our simulations. The result appears robust despite the large value of significant wave height, Hs/za ; 0.64. Examination of the results over the interval 0 , z , za shows that the total stress at the surface is slightly larger (more negative) than at za. A small vertical flux divergence for both u and y components, about 8% of u2 *,a, is required to balance the Coriolis forces and the large-scale pressure gradients that appear in (8c). A noticeable decrease with height in the u momentum fluxes starts at about z/zi ' 0.1. The near-constant-flux layer observed in Fig. 10 is not imposed but is a consequence of the fine mesh and resolved dynamics in the LES. Since the present SGS model is an adaption of a standard TKE model (Moeng and Wyngaard 1988), we expect other comparable SGS models will generate similar solutions on a fine mesh.

Figures 10 and 11 also illustrate important wave effects. The partitioning of the total momentum flux between tf, pf, and sf is strongly height and waveage dependent. For example, in the simulation with Cp/Ua 5 4.2, wave effects are readily apparent in the turbulent and pressure flux contributions well above za extending as high as 100m or more (also see Figs. 17 and 18); on average, they tend to counterbalance each other for z . za. The shifting flux balance between tf, pf, and sf below za is a consequence of the scale content of the wave spectrum and the magnitude of the surface winds. For the simulation with Cp/Ua 5 4.2, the longest wavelength modes in the wave spectrum propagate rapidly compared to the surface winds and at z 5 0 generate net upward momentum flux through a positive pressure-wave slope correlation. This positive form stress is counterbalanced by a large negative SGS contribution resulting in a surface stress that is nearly equal to its value at za. Note that the resolved turbulent flux equals zero at z50 because of no flow through the water surface. The implication for this large wave-age simulation is that small-scale waves, smaller than are resolved by the LES mesh, slow the winds at z 5 0. As we decrease the wave age, by increasing the geostrophic winds, the momentum flux partitioning smoothly shifts between tf, pf, and sf. The large-scale waves support small form stress as they are propagating at a speed close to the surface winds, and then the intermediate- and smaller-scale wave modes play an important role in supporting the wind stress through a negative pressure- wave slope correlation (form drag). The overall level of the resolved form stress compared to the total wind stress in the LES hints that the high-wavenumber tail of the wave spectrum plays an important role in generating drag at wind-wave equilibrium conditions. The spectral content of the form stress is given by the cospectrum between surface pressure p* and wave slope zx/J shown in Fig. 12. There is significant wave-age variability with a local peak and sign changes at low wavenumbers, kx , 0.4m21. Högström et al. (2009) also suggests that a decomposition of the total drag will include both positive and negative contributions because of swell. However, at high wavenumbers, the small scales, say l,15 m, add a negative contribution to the overall form stress out to the largest resolved wavenumber in the LES for all wave ages considered (see the bottom panel of Fig. 12). Makin et al. (1995, see their Fig. 4), using a wave growth model, estimate that scales between 0.1 to 10m carry a large fraction of the form drag.

Figure 13 shows the profile of the total heat flux over the complete MABL. Key points to notice are the dramatic increase in entrainment flux at the MABL top with increasing wind speeds, which is a direct consequence of the large entrainment events observed in the flow visualization of Fig. 6. The (negative) normalized entrainment flux increases from 0.12 to 1.09 as the geostrophic winds increase from 5 to 25ms21, with the height of the minimum entrainment flux descending with increasing wind speed; the width of the entrainment zone, loosely defined as region where the scalar flux is negative, also expands. The enhancement of entrainment flux at near-neutral conditions is significant and we find similar and larger increases in entrainment flux in ocean boundary layers driven by hurricane winds and waves (Sullivan et al. 2012). Thus, for the high-wind MABLs examined here, scalar transport is likely influenced by entrainment dynamics. Near and below za the total scalar flux is nearly constant for varying wave age, adding support to the assumption of a constant-flux layer above surface waves.

The resolved turbulence variances s2 u,y,w shown in Fig. 14 further illustrate wind-wave couplings that vary with wave age: variance is defined as s2f 5h[ f (j, h)2hf i]2i. Notice in the simulation with large wave age, so-called ''wave pumping'' (McWilliams et al. 1999; Sullivan et al. 2008) by the long wavelength modes in the wave spectrum is vigorous and greatly enhances the variances generating maxima near the water surface. The surfaceenhanced variance decays slowly with increasing z extending upward to z/zi .0.1. We speculate that nearsurface wave pumping might lead to an enhancement of the horizontal u- and y-variance components through pressure-strain correlations that are opposite to the usual energy redistribution in a flat-wall boundary layer [see p. 105 in Wyngaard (2010)]. At wind-wave equilibrium, the background shear-generated turbulence overpowers the small-scale wave-induced motions and then the shape and magnitude of the variance profiles are qualitatively similar to those obtained over stationary waves. Also the w variance in the simulation with Cp/Ua 5 4.2 reveals a local maximum slightly below z/zi;0.5 as expected for a boundary layer forced by convection.

e. Wind-wave correlations Wind-wave studies seek to identify wave-induced (correlated) motions in the turbulent surface-layer winds as this provides insight into wave growth (e.g., Belcher and Hunt 1998; Hristov et al. 1998, 2003; Sajjadi et al. 2014).Aqualitative indication of the instantaneous correlation between vertical velocity and the wave field for varying wave age is shown in the snapshots in Figs. 15 and 16; the correlations are quantified in Figs. 17 and 18. The images reflect the atmospheric response for varying wind speeds at a fixed height when the wave spectrum is held constant. At the reference height za, the amplitude and organization of the wave signal in vertical velocity is clearly wave-age dependent. Despite the randomness of the wind and wave fields, one can readily identify the imprint of the wave field on vertical velocity and notice its decay with decreasing wave age. This is expected based on the flow statistics discussed previously and the visualization presented in section 4b. Close to windwave equilibrium (Cp/Ua 5 1.1), there is only a slight hint of the wave signature in the vertical velocity field. As Cp/Ua increases past equilibrium, the wave signature gradually emerges and eventually dominates. Past simulations and observations in idealized calculations hint at this dependence and the present calculations with a realistic wave spectrum confirm this wave-age dependence.

Figures 17 and 18 show vertical profiles of the windwave correlations in our computations. The correlation coefficients are defined as ... (32) where sf denotes the standard deviation of variable f. The correlations illustrate intriguing wave-age dependencies. In large wave-age (nonequilibrium) conditions, the large wave modes propagate quickly compared to the surface winds and their imprint on the atmosphere is best observed in the correlation between vertical velocity and the local time variation of the wave height (i.e., hw?h/?ti). In other words, vertical velocity is nearly 908 out of phase with respect to the wave surface. Note that a finite correlation between w and ?h/?t extends well above za and is a clear indication of wave pumping on the atmospheric winds by the large-amplitude, fast-movingwavemodes. This strong correlation disappears near wind-wave equilibrium and subsequently vertical velocity becomes more weakly correlated with wave height in a shallow layer well below za. The normalized correlation between u and (h, ?h/?t) is opposite to that of the vertical velocity, as expected. The highest correlation occurs in the swell-dominated regime with wave height. The correlation between u and the time change of wave height is complex and confined to the layer below za and may be due to critical-layer effects as reported by Hristov et al. (1998, 2003) and Grare et al. (2013).Hence, detailedmeasurements in a wave-following coordinate system very near thewater surface are required to extract wave-induced momentum flux from observational data under wind-wave equilibrium conditions.

5. Summary and conclusions A large-eddy simulation (LES) model that adopts a Boussinesq flow model for the marine atmospheric boundary layer (MABL) coupled to a broadband timedependent surface gravity wave field is described. The algorithmuses a coordinate transformation from physical to computational space that is wave following and accounts for the vertical motion of themesh (or grid speed) in physical space: the time dependence of the grid and the grid speeds appear explicitly in themomentumand scalar time tendency and advection terms. The LES model is used to examine the wind-wave dynamics in a weakly unstable MABL for a range of geostrophic wind speeds Ug 5 5, 10, 20, and 25ms21 with surface heat flux Q* 5 0.01Kms21. The broadband moving wave field is taken to be a wind-wave equilibrium spectrum with the significant wave heightHs56.4m. The wave age Cp/Ua5[1.1, 4.2] varies from near wind-wave equilibrium to a lowwind swell-dominated regime.

The high-resolution simulations show that nearsurface wave-impacted turbulence smoothly blends into the bulk of the MABL over a vertical distance that depends on wave age. For low winds and strong swell with wave age (Cp/Ua 5 4.2), wave signals can be found over an extended vertical distance. In this situation, wind speeds normalized by friction velocity computed at the standard reference height (za510 m) are faster than their counterparts over a rough surface with no resolved waves. As the winds and waves approach equilibrium and the background turbulence increases, wave-induced signals are mainly generated by small-scale waves and their wave signatures are confined close to the water surface below za. Wind profiles from Monin-Obukhov (MO) similarity formulas deviate from the LES wind profiles in the situation with large wave age. At windwave equilibrium, the agreement between LES andMO is good over a vertical range extending up to z/zi ; 0.2. However, even at wind-wave equilibrium conditions, flow over moving waves differs from its counterpart over an identical but stationary wave spectrum. In the situation of stationary bumps, the wavy surface has a larger roughness parameter and the vertical height where the winds blend into a log-linear profile is higher.

A nearly-constant-flux layer for streamwise momentum extends up to za; the momentum flux budget computed in wave-following coordinates includes turbulent flux, pressure-wave slope correlations, and subgrid-scale contributions. A similar budget for spanwise momentum flux decreases by about 8%over 0 , z , za. Belowza the momentum flux balance smoothly shifts between the turbulent flux, pressure-wave slope, and SGS contributions as the water surface is approached. At wind-wave equilibrium, the wave-induced formstress is negative and primarily confined close to the water surface while, at large wave age, wave-driven winds are produced by a positive form stress that extends well above za. The cospectrum between surface pressure and wave slope, essentially the spectral content of the form stress, switches sign between positive and negative values at low wavenumbers (kx , 0.4m21) depending on wave age. At higher wavenumbers, the form stress is negative and small scales (l , 15m) contribute to the overall surface drag. This suggests that higher-resolution LESs are needed to evaluate the viscous and pressure contributions to surface drag at wind-wave equilibrium. The resolved velocity variances also exhibit wave-age dependence. Higher variances, especially for vertical velocity, are produced by wave pumping at large wave age.

Flow visualization and wind-wave correlations further illustrate the dependence on wave age. Vertical (w) and streamwise (u) velocity are best correlated with ?h/?t and h, respectively; at large wave age, the correlations approach (1, 21) at the water surface. At wind-wave equilibrium, w and u are only weakly correlated with wave height at za. For the combination of weak surface heating and winds considered, depth-filling shear-convective rolls develop in theMABL and the entrainment flux at the top of the boundary layer depends strongly on wind speed. Large entrainment events extend from the inversion down to the middle of the MABL.

Acknowledgments. PPS and JCM were supported by the Physical Oceanography Program through the Office of Naval Research (ONR) Departmental Research Initiative, High-Resolution Air-Sea Interaction (HiRes). PPS and EGP are supported by the Department of Energy through the Wind and Water Power Technology Office (DE-EE0005373) and also by the National Science Foundation through the National Center for Atmospheric Research. We thank the many HiRes investigators, Tetsu Hara, and the three anonymous reviewers who provided insights that helped to improve this work. This research benefited from computer resources provided by the NCAR Strategic Capability program and the Department of Defense High-Performance Computing Modernization Program. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

1 Thomas and Lombard (1979, p. 1032) develop the geometric conservation law using an integral formulation but note that it can equally be derived from the differential form of the governing equations by setting r 5 1 and ui 5 0, which is the procedure used here.

REFERENCES Anderson, D. A., J. C. Tannehill, and R. H. Pletcher, 1984: Computational Fluid Mechanics and Heat Transfer. McGraw-Hill, 599 pp.

Banner, M. L., 1990: The influence of wave breaking on the surface pressure distribution in wind-wave interaction. J. Fluid Mech., 211, 463-495, doi:10.1017/S0022112090001653.

_____, and W. K. Melville, 1976: On the separation of airflow over water waves. J. Fluid Mech., 77, 825-842, doi:10.1017/ S0022112076002905.

Belcher, S. E., and J. C. R. Hunt, 1998: Turbulent flow over hills and waves. Annu. Rev. Fluid Mech., 30, 507-538, doi:10.1146/ annurev.fluid.30.1.507.

Bell, M. M., M. T. Montgomery, and K. A. Emanuel, 2012: Air-sea enthalpy and momentum exchange at major hurricane wind speeds observed during CBLAST. J. Atmos. Sci., 69, 3197- 3222, doi:10.1175/JAS-D-11-0276.1.

Chalikov, D. V., 1978: The numerical simulation of wind-wave interaction. J. Fluid Mech., 87, 561-582, doi:10.1017/ S0022112078001767.

_____, 1998: Interactive modeling of surface waves and boundary layer. Ocean Wave Measurement and Analysis: Proceedings of the Third International Symposium WAVES 97, B. L. Edge, J. M. Hemsley, and Y. Goda, Eds., American Society of Civil Engineers, 1525-1539.

Chen, S. S., W. Zhao, M. A. Donelan, and H. L. Tolman, 2013: Directional wind-wave coupling in fully coupled atmosphere- wave-ocean models: Results from CBLAST-Hurricane. J. Atmos. Sci., 70, 3198-3215, doi:10.1175/JAS-D-12-0157.1.

Deardorff, J. W., 1972: Three-dimensional numerical modeling of the planetary boundary layer. Workshop on Micrometeorology, D. A. Haugen, Ed., Amer. Meteor. Soc., 271-311.

Demirdzic, I., and M. Peric, 1990: Finite volume method for prediction of fluid flow in arbitrarily shaped domains with moving boundaries. Int. J. Numer. Methods Fluids, 10, 771-790, doi:10.1002/fld.1650100705.

Donelan, M. A., 1998: Air-water exchange processes. Physical Processes in Lakes and Oceans, J. Imberger, Ed., Coastal and Estuarine Studies, Vol. 54, Amer. Geophys. Union, 19-36.

_____, J. Hamilton, and W. H. Hui, 1985: Directional spectra of wind-generated waves. Philos. Trans. Roy. Soc. London, 315A, 509-562, doi:10.1098/rsta.1985.0054.

_____, B. K. Haus, N. Reul, W. J. Plant, M. Stiassnie, H. C. Graber, O. B. Brown, and E. S. Saltzman, 2004: On the limiting aerodynamic roughness of the ocean in very strong winds. Geophys. Res. Lett., 31, L18306, doi:10.1029/2004GL019460.

Druzhinin, O., Y. Troitskaya, and S. Zilitinkevich, 2012: Direct numerical simulation of a turbulent wind over a wavy water surface. J. Geophys. Res., 117, C00J05, doi:10.1029/2011JC007789.

Edson, J., and Coauthors, 2007: The coupled boundary layers and air-sea transfer experiment in low winds. Bull. Amer. Meteor. Soc., 88, 341-356, doi:10.1175/BAMS-88-3-341.

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, doi:10.1175/1520-0442(2003)016,0571:BPOASF.2.0.CO;2.

Fan, Y., S.-J. Lin, I. M. Held, Z. Yu, and H. L. Tolman, 2012: Global ocean surface wave simulation using a coupled atmosphere-wave model. J. Climate, 25, 6233-6252, doi:10.1175/ JCLI-D-11-00621.1.

Fedorovich, E., and Coauthors, 2004: Entrainment into sheared convective boundary layers as predicted by different large eddy simulation codes. 16th Symp. on Boundary Layers and Turbulence, Portland, ME, Amer. Meteor. Soc., P4.7. [Available online at https://ams.confex.com/ams/BLTAIRSE/ techprogram/paper_78656.htm.] Gent, P. R., and P. A. Taylor, 1976: A numerical model of the air flow above water waves. J. Fluid Mech., 77, 105-128, doi:10.1017/S0022112076001158.

Ghosal, S., and P. Moin, 1995: The basic equations for the large eddy simulation of turbulent flows in complex geometry. J. Comput. Phys., 118, 24-37, doi:10.1006/jcph.1995.1077.

Goldstein, H., 1950: Classical Mechanics. Addison-Wesley Publishing Company, 399 pp.

Grare, L., L. Lenain, and W. Melville, 2013: Wave-coherent airflow and critical layers over ocean waves. J. Phys. Oceanogr., 43, 2156-2172, doi:10.1175/JPO-D-13-056.1.

Hanley, K. E., and S. E. Belcher, 2008: Wave-driven wind jets in the marine atmospheric boundary layer. J. Atmos. Sci., 65, 2646- 2660, doi:10.1175/2007JAS2562.1.

_____, _____, and P. P. Sullivan, 2010: A global climatology of wind- wave interaction. J. Phys. Oceanogr., 40, 1263-1282, doi:10.1175/2010JPO4377.1.

Henn, D. S., and R. I. Sykes, 1999: Large-eddy simulations of flow over wavy surfaces. J. Fluid Mech., 383, 75-112, doi:10.1017/ S0022112098003723.

Högström, U., A. Smedman, E. Sahlée, W. M. Drennan, K. K. Kahma, H. Pettersson, and F. Zhang, 2009: The atmospheric boundary layer during swell: A field study and interpretation of the turbulent kinetic energy budget for high wave ages. J. Atmos. Sci., 66, 2764-2779, doi:10.1175/2009JAS2973.1.

_____, A. Rutgersson, E. Sahlée, A. Smedman, T. Hristov, W. M. Drennan, and K. K. Kahma, 2013: Air-sea interaction features in the Baltic Sea and at a Pacific trade-wind site: An intercomparison study. Bound.-Layer Meteor., 147, 139-163, doi:10.1007/s10546-012-9776-8.

Hristov, T., C. Friehe, and S. Miller, 1998: Wave-coherent fields in air flow over ocean waves: Identification of cooperative behavior buried in turbulence. Phys. Rev. Lett., 81, 5245-5248, doi:10.1103/PhysRevLett.81.5245.

_____, S.D.Miller, and C. Friehe, 2003: Dynamical coupling of wind and ocean waves through wave-induced air flow. Nature, 422, 55-58, doi:10.1038/nature01382.

Komen, G. J., L. Cavaleri, M. Donelan, K. Hasselmann, S. Hasselmann, and P. A. E. M. Janssen, 1994: Dynamics and Modelling ofOceanWaves. Cambridge University Press, 532 pp.

Large, W. G., J. C. McWilliams, and S. C. Doney, 1994: Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. Rev. Geophys., 32, 363-403, doi:10.1029/94RG01872.

Le, H., and P. Moin, 1991: An improvement of fractional step methods for the incompressible Navier-Stokes equations. J. Comput. Phys., 92, 369-379, doi:10.1016/0021-9991(91)90215-7.

Li, P. Y., 1995: A numerical study on energy transfer between turbulent air flow and finite amplitude water waves. Ph.D. thesis, York University, North York, Ontario, Canada, 182 pp.

Lighthill, M. J., 1962: Physical interpretation of the mathematical theory of wave generation by wind. J. Fluid Mech., 14, 385- 398, doi:10.1017/S0022112062001305.

Makin, V. K., V. N. Kudryavtsev, and C. Mastenbroek, 1995: Drag of the sea surface. Bound.-Layer Meteor., 73, 159-182, doi:10.1007/BF00708935.

McWilliams, J. C., C.-H. Moeng, and P. P. Sullivan, 1999: Turbulent fluxes and coherent structures in marine boundary layers: Investigations by large-eddy simulation. Air-Sea Exchange: Physics, Chemistry, Dynamics, and Statistics, G. Geernaert, Ed., Kluwer, 507-538.

Miles, J. W., 1957: On the generation of surface waves by shear flow, Part I. J. Fluid Mech., 3, 185-204, doi:10.1017/ S0022112057000567.

Moeng, C.-H., 1984: A large-eddy simulation model for the study of planetary boundary-layer turbulence. J. Atmos. Sci., 41, 2052-2062, doi:10.1175/1520-0469(1984)041,2052: ALESMF.2.0.CO;2.

-, and J. C. Wyngaard, 1988: Spectral analysis of large-eddy simulations of the convective boundary layer. J. Atmos. Sci., 45, 3573-3587, doi:10.1175/1520-0469(1988)045,3573: SAOLES.2.0.CO;2.

-, and P. P. Sullivan, 1994: A comparison of shear and buoyancy driven planetary-boundary-layer flows. J. Atmos. Sci., 51, 999-1022, doi:10.1175/1520-0469(1994)051,0999: ACOSAB.2.0.CO;2.

Nakayama, A., H. Noda, and K. Maeda, 2004: Similarity of instantaneous and filtered velocity fields in the near wall region of zero-pressure gradient boundary layer. Fluid Dyn. Res., 35, 299-321, doi:10.1016/j.fluiddyn.2004.07.002.

-,K. Sakio, Y. Kitano, and S. Yokojima, 2009: Direct numerical simulation and filtering of turbulent flows over model rough surfaces. Mem. Grad. School Eng. Kobe Univ., 1, 9-41.

Nilsson, E. O., A. Rutgersson, A. Smedman, and P. Sullivan, 2012: Convective boundary layer structure in the presence of windfollowing swell. Quart. J. Roy. Meteor. Soc., 138, 1476-1489, doi:10.1002/qj.1898.

Phillips, O. M., 1977: Dynamics of the Upper Ocean. Cambridge University Press, 336 pp.

Richter, D. H., and P. P. Sullivan, 2013: Sea surface drag and the role of spray. Geophys. Res. Lett., 40, 656-660, doi:10.1002/ grl.50163.

Romero, L., and W. K. Melville, 2010: Airborne observations of fetch-limited waves in the Gulf of Tehuantepec. J. Phys. Oceanogr., 40, 441-465, doi:10.1175/2009JPO4127.1.

Sajjadi, S. G., J. C. R. Hunt, and F. Drullion, 2014: Asymptotic multi-layer analysis of wind over unsteady monochromatic surface waves. J. Eng. Math., 84, 73-85, doi:10.1007/ s10665-013-9663-4.

Spalart, P. R., R. D. Moser, and M. M. Rogers, 1991: Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions. J. Comput. Phys., 96, 297-324, doi:10.1016/0021-9991(91)90238-G.

Sullivan, P. P., and J. C. McWilliams, 2010: Dynamics of winds and currents coupled to surface waves. Annu. Rev. Fluid Mech., 42, 19-42, doi:10.1146/annurev-fluid-121108-145541.

-, and E. G. Patton, 2011: The effect of mesh resolution on convective boundary-layer statistics and structures generated by large-eddy simulation. J. Atmos. Sci., 68, 2395-2415, doi:10.1175/JAS-D-10-05010.1.

-, J. C. McWilliams, and C.-H. Moeng, 1996: A grid nesting method for large-eddy simulation of planetary boundary layer flows. Bound.-Layer Meteor., 80, 167-202, doi:10.1007/ BF00119016.

-, C.-H. Moeng, B. Stevens, D. H. Lenschow, and S. D. Mayor, 1998: Structure of the entrainment zone capping the convective atmospheric boundary layer. J. Atmos. Sci., 55, 3042-3064, doi:10.1175/1520-0469(1998)055,3042:SOTEZC.2.0.CO;2.

-, J. C. McWilliams, and C.-H. Moeng, 2000: Simulation of turbulent flow over idealized water waves. J. Fluid Mech., 404, 47-85, doi:10.1017/S0022112099006965.

-, T. W. Horst, D. H. Lenschow, C.-H. Moeng, and J. C. Weil, 2003: Structure of subfilter-scale fluxes in the atmospheric surface layer with application to large-eddy simulation modeling. J. Fluid Mech., 482, 101-139, doi:10.1017/ S0022112003004099.

-, J. B. Edson, T. Hristov, and J. C. McWilliams, 2008: Largeeddy simulations and observations of atmospheric marine boundary layers above nonequilibrium surface waves. J. Atmos. Sci., 65, 1225-1245, doi:10.1175/2007JAS2427.1.

-, L. Romero, J. C. McWilliams, and W. K. Melville, 2012: Transient evolution of Langmuir turbulence in ocean boundary layers driven by hurricane winds and waves. J. Phys. Oceanogr., 42, 1959-1980, doi:10.1175/JPO-D-12-025.1.

Suzuki, N., T. Hara, and P. P. Sullivan, 2013: Impact of breaking wave form drag on near-surface turbulence and drag coefficient over young seas at high winds. J. Phys. Oceanogr., 43, 324-343, doi:10.1175/JPO-D-12-0127.1.

Thomas, P. D., and C. K. Lombard, 1979: Geometric conservation law and its application to flow computations on moving grids. AIAA J., 17, 1030-1037, doi:10.2514/3.61273.

Veron, F., G. Saxena, and S. K. Misra, 2007: Measurements of the viscous tangential stress in the airflow above wind waves. Geophys. Res. Lett., 34, L19603, doi:10.1029/2007GL031242.

Wyngaard, J. C., 2004: Toward numerical modeling in the Terra Incognita. J. Atmos. Sci., 61, 1816-1826, doi:10.1175/ 1520-0469(2004)061,1816:TNMITT.2.0.CO;2.

_____, 2010: Turbulence in the Atmosphere. Cambridge University Press, 393 pp.

Yang, D., and L. Shen, 2010: Direct-simulation-based study of turbulent flow over various waving boundaries. J. Fluid Mech., 650, 131-180, doi:10.1017/S0022112009993557.

_____, C. Meneveau, and L. Shen, 2013: Dynamic modelling of seasurface roughness for large-eddy simulation of wind over ocean wavefield. J. Fluid Mech., 726, 62-99, doi:10.1017/ jfm.2013.215.

Zang, Y., R. L. Street, and J. R. Koseff, 1994:Anon-staggered grid, fractional step method for time-dependent incompressible Navier-Stokes equations in curvilinear coordinates. J. Comput. Phys., 114, 18-33, doi:10.1006/jcph.1994.1146.

Zheng, Z., and L. Petzold, 2006: Runge-Kutta-Chebyshev projection method. J. Comput. Phys., 219, 976-991, doi:10.1016/ j.jcp.2006.07.005.

PETER P. SULLIVAN National Center for Atmospheric Research,* Boulder, Colorado JAMES C. MCWILLIAMS Department of Atmospheric and Oceanic Science, University of California, Los Angeles, Los Angeles, California EDWARD G. PATTON National Center for Atmospheric Research,* Boulder, Colorado (Manuscript received 21 April 2014, in final form 1 July 2014) * The National Center for Atmospheric Research is sponsored by the National Science Foundation.

Corresponding author address: Peter P. Sullivan, Mesoscale and Microscale Meteorology Division, P.O. Box 3000, National Center for Atmospheric Research, Boulder, CO 80307-3000.

E-mail: [email protected] (ProQuest: Appendix omitted.) (c) 2014 American Meteorological Society

[ Back To TMCnet.com's Homepage ]