Wave-Modiﬁed Ekman Current Solutions for the Time-Dependent Vertical Eddy Viscosity

: This study aimed to highlight a general lack of clarity regarding the scale of the temporal averaging implicit in Ekman-type models. Under the assumption of time and depth-dependent eddy viscosity, we present an analytical Fourier series solution for a wave-modiﬁed Ekman model. The depth dependence of eddy viscosity is based on the K-Proﬁle Parameterization (KPP) scheme. The solution reproduces major characteristics of diurnal variation in ocean velocity and shear. Results show that the time variability in eddy viscosity leads to an enhanced mean current near-surface and a decrease in the effective eddy viscosity, which ﬁnally results in an intensiﬁed near-surface shear and wakes a low-level jet ﬂow. Rectiﬁcation values are dominated by the strength of diurnal mixing, and partly due to the nonlinear depth dependence of the eddy viscosity.


Introduction
Accuracy of the vertical structure of near-surface flow is of prime importance for modeling the horizontal transport of pollutants, algae, chlorophyll, and so on. Though there are some discrepancies between the observation and classical Ekman theory, the latter is still the foundation of our understanding of ocean surface currents. The turbulent mixing in the ocean boundary layer is known to vary with time, and many affected physical processes, such as surface gravity waves, air entrainment, solar heating, density stratification, and Langmuir circulations are not fully understood yet [1][2][3][4].
Diurnal heat flux between the air-sea boundary layer results in a daily cycle of temperature, stratification, and mixing in the upper ocean, which probably leads to a diurnal cycle in eddy viscosity. Recent observations declare the importance of diurnal variations of stratification and eddy viscosity [1,2,5,6]. Shrira and Almelah [7] suggested that the Ekman theory should be extended to take into account time-and depth-dependent eddy viscosity, which is expected to be a better reflection of reality. McWilliams et al. [8] used Large Eddy Simulation (LES) to prove that the effective eddy viscosity should not be horizontally uniform or isotropic, but anisotropic vectors.
Observations proved the effects of surface gravity waves cannot be neglected in the ocean surface current. Following the idea that the wave-averaged Stokes drift should also be taken into account in the Ekman model by Huang [9], the wave-modified Ekman model was analyzed theoretically [10][11][12][13][14][15] and numerically by wave-ocean coupled models [16,17] or LES [3,18,19].
Most previous studies about the ocean surface boundary layer are based on the assumptions of independence of time of both eddy viscosity and the Stokes drift. What happens without this assumption is that one of the issues need to be interpretation. This work was inspired by the study of Wenegrat and McPhaden [6], which applied the approaches used extensively in the dynamics of low-level jets in the atmospheric boundary layer [20][21][22] to oceanography. As pointed out in Wenegrat and McPhaden [6], understanding the effect of the time variability in eddy viscosity is critical to proper use of timeaveraged observations, since observations are often averaged in different time-scales in order to remove noise contamination. Under the assumption of time-and depth-dependent eddy viscosity, we propose a Fourier series solution for a non-stationary Ekman problem, discuss the qualitative characteristics of this solution, and evaluate the rectification of diurnal eddy viscosity.
The paper is organized as follows. In Section 2, we provide the mathematical model, then a Fourier series solution for the time-and depth-dependent viscosity is derived in Section 3. Section 4 is mainly about the discussion of the solution: the sensitivity to Fourier modes is firstly tested, then a scenario of a diurnal cycle of KPP eddy viscosity in a wavemodified Ekman model is presented; and finally, the rectification effect is conducted. Our findings and conclusions are summarized in the last section.

Basic Equations and Boundary Conditions
When the effects of random surface waves are considered, with zero horizontal pressure gradient and advection and no momentum flux through the lower boundary layer, the wind-driven horizontal current satisfies the following modified momentum equation [15,17,23,24] if we assume that the wave, velocity, and turbulent properties are uniform horizontally: where U = U(z, t) = u(z, t) + iv(z, t) is the complex horizontal velocity in the x − y plane, and the horizontal coordinate axes are fixed on the still water level with z = 0. The schematic diagram of the coordinate system is shown in Figure 1, in which we assumed the wave direction is the same as the wind stress. However, in nature, the wave direction does not need to be the same as the wind, and there always exists a wind-wave misalignment. i = √ −1, f cor is the vertical component of the Coriolis frequency and t is time, A v is the vertical eddy viscosity, U s = u s + iv s is the complex Stokes drift, and T wds is the wave-induced momentum transfer from waves to the mean flow due to the dissipation of wave energy. The velocity U discussed here is the quasi-Eulerian current, which is equal to the Lagrangian mean current minus the Stokes drift and can be understood as the Eulerian-mean current, as stated by Jenkins [23,25] and Perrie et al. [16]. The vertical eddy viscosity is generally a function of both time and space. To better explain the effects of the near-surface ocean diurnal cycle on horizontal currents in the ocean surface boundary layer, a time-varying turbulent viscosity is considered. Following Zhang and Tan [17] and Wenegrat and McPhaden [6], assuming the vertical eddy viscosity varies with both time and depth, and the dependence on time and depth is non-related, then it can be approximately written as which ensures the known eddy viscosity can be separable in time and space. The time dependence takes a particular form where ω = 2π/86,400 stands for the diurnal frequency, which also equals to the angular velocity of the Earth, and δ ∈ [0, 1) can be adjusted to determine the strength of the periodic cycle of mixing. This assumption was commonly used in the study of the dynamics of lowlevel jets in the atmospheric boundary layer [20,22,26,27]. Applied to the oceanographic problem by Wenegrat and McPhaden [6], based on the observational data at the tropical Atlantic from 13 October 2008 to 6 January 2009, which obtained the diurnal cycle of nearsurface eddy viscosity [5] and suggested a sinusoidal time-dependence is a reasonable first approximation, as shown in Figure 1 of Wenegrat and MaPhaden [6]. This assumption was accepted by other authors and used in diurnal cycling of submesoscale dynamics [28][29][30].
For the depth dependence of eddy viscosity A(z), it can be parameterized as the simplified K-Profile Parameterization (KPP) scheme [1,8,15,31], which can be formulated as the following: In which u * is the oceanic friction velocity, h b is the boundary layer depth, and c 1 and c 2 are two constants. The normalized depth σ increases from 0 at the surface to 1 at the bottom of the boundary layer. The values of c 1 and c 2 are taken as c 1 = κ = 0.4 to be consistent with a wall-bounded similarity layer (i.e., log layer), where A(z) → −κu * z as z → 0, and c 2 is taken as 2.0, which is obtained by fitting the viscosity profiles of the LES presented by McWilliams et al. [8]. As shown by Song and Xu [15], the wave-modified Ekman current solutions calculated with the choice of c 1 = 0.4 and c 2 = 2.0 exhibit good agreement with those well-known published observational data of the Ekman layer, and the solutions are not sensitive to these two constant parameters. By the way, in order to avoid zero eddy viscosity, the vertical coordinate always starts from z = −1 m to the rounding of boundary layer depth.
The surface boundary condition for the wave-modified Ekman current is where ρ w is the water density, τ in is the reduction of wind stress due to wave generations, and τ a is the complex wind stress, computed from surface wind field U 10 at 10 m height, where ρ a is the air density,θ * is the unit vector in that direction, and C d is the air-sea drag coefficient, which is related to U 10 by the relation [32] C d = (0.8 + 0.065|U 10 |) × 10 −3 .
The lower boundary condition is taken as a no-slip bottom boundary Generally, the lower boundary is assumed to lie well below the surface boundary layer. Here, we assume that H = −h b . The influence of this assumption has been examined with different c 2 ; it will not change the pattern of the results shown below, and it generally satisfies the bottom boundary condition in our model.
From the above momentum, for Equation (1), along with boundary condition (5), three wave-related terms were included in the Ekman model, namely, the Coriolis-Stokes Force (i f cor U s , here after CSF), the reduction of wind stress due to wave generations (τ in ), and momentum transfer from waves to the mean flow due to wave dissipation T wds , in which τ in and T wds can be estimated by the source terms from a directional spectral wave prediction model, which act to transfer momentum from wind to wave and wave to ocean, as follows [13,17,23]: where f is the frequency of the wave, k is the moduli of the horizontal wavenumber K and K = k x + ik y = k cos θ + ik sin θ, their relationship is given by the deep water dispersion relation ω 2 = (2π f ) 2 = gk, and θ is the angle between the wave vector and the x-axis. S in ( f , θ) is the source term of wind input energy to waves, and S ds ( f , θ) is the source term of wave energy loss due to wave dissipation; details of the related mechanisms and parameterization schemes are described in third-generation WAM-type models [33,34]. Stokes drift U s can be expressed as in Kenyon [35] where E( f , θ) is the directional frequency spectrum of surface waves.

Solutions
Transforming where Using (2) and (3), and transforming the time coordinate, such that where Substituting (4) into (14) and (15), the wave-modified equation of (1) can be written as where Then, the boundary conditions (5) and (8) become where and Following Song and Xu [15], a complex Fourier series solution can be founded by expanding G(x, ζ) and forcing variables D(x, ζ), τ(ζ) into Fourier series where ν n = n/T ζ is the discrete frequency, T ζ is the period of ζ, and G n (x), D n (x), and τ n are the complex Fourier coefficients of G(x, ζ), D(x, ζ), and τ(ζ) for the discrete frequency ν n , which are independent with time and can be determined by Here, υ n = n/T and T is the period of the diurnal cycle. Substituting (22)-(24) into (16)- (21), the nth component of the nth momentum equation of (16) can be written as with the boundary conditions Note that (28)- (30) are the same as the results of Song and Xu [15]. Hence, the general solution of Equation (28) can be given by where Here, and F(a, b; c; x) is the hypergeometric function in (39) and (40). According to the above derivation, the wave-modified Ekman current for the timedependent vertical eddy viscosity formulated by the KPP scheme (4) has the Fourier series solution as follows:

Discussion of the Solution
Since we focus on the specific diurnal cycle of eddy viscosity here, we do not consider the transient wind, but simply assume steady wind forcing. Along with steady wind forcing, we suppose that the wave is fully developed and that the mean wave propagates in the same direction as the wind stress. We use the commonly monochromatic approximation for Stokes drift profiles U s (z) = U s0 exp(z/h s ), since it had been effectively used in theoretical illustration and idealized numerical modeling [5,12,36,37], where U s0 is the surface Stokes drift velocity and h s is the Stokes depth scale. τ in and T wds were ignored here, as their effects were much smaller than CSF, for there is a counterbalance between the wind-input-wave stress and wave-to-ocean dissipation for a fully developed wind wave [13,24].
The validation of the KPP scheme in the Ekman model has been discussed and compared with some numerical simulations and observations [1,2,38]; therefore, we believe it can be used to provide a reasonable result for the problem of rectification in an ocean surface boundary layer under a time-varying eddy viscosity. For the eddy viscosity in the KPP scheme, (4) has a dependence on the latitude, and is expected to give a difference rectification characteristic.

Sensitivity to the Fourier Modes n
As the solution given by (42) is related to Fourier modes n, which is closely related to the accuracy of the results, it is of critical importance to determine a reasonable order of n in a different strength of the periodic cycle of mixing. Because an order of n that is too high would largely increase the computational burden, while an order that is too low could bring about irregular errors for the surface current, an economic cut-off order of n needs to be examined and selected during the calculation. As shown in Figure 2, under the diurnal cycle of eddy viscosity, for the same periodic strength as δ = 0.7, the surface zonal velocity u 0 displays different degrees of fluctuation at the first few lower Fourier modes, when n = ±5, ∂u 0 /∂t fluctuates over time; while smooth u 0 is shown in the last few high modes, ∂u 0 /∂t continuously changes and shows a daily trend. With the increase of δ, the minimum of the Fourier modes n resulted in a smooth time-varying velocity getting higher, and for δ = 0.9, the order of Fourier modes needs to be larger than 55, while for δ = 0.3, n = ±5 is sufficient to give a rather stable and accurate result. Therefore, the following contents are based on reasonable and effective Fourier series modes. The latitude will also influence the flatness of the surface velocity varying with time; as shown in formulation (42), both the eddy viscosity in the KPP scheme and CSF are related to the latitude. However, it will not significantly influence the critical cut-off of n for the same δ. As shown in Figure 3, a slight increase of the order of n will be sufficient to ensure validity for all the latitudes.
It is worth mentioning that the Fourier series solution will also encounter a similar problem, as discussed in Appendix B of Wenegrat and McPhaden [6]. When f cor + 2πν n = 0, which happens when n = ±1, 2 at the latitude 30°and 90°, this will lead to a zero denominator for A n1 and A n2 in (34) and (35). However, as shown in Figure 3b, when the latitude equaled 30°, the Fourier series solution still performed well. It again confirmed that the contribution of these specific modes to the total solution is rather small, as interpreted in Wenegrat and McPhaden [6].

Qualitative Characteristics
Under middle latitude and forced by a constant zonal wind stress of U 10 = 10 m/s, Figure 4 gives an example of the time-averaged velocity profile with diurnal eddy viscosity (δ = 0.6) compared with no time-varying eddy viscosity (δ = 0). The dotted lines in Figure 4 stand for δ = 0 without considering CSF, and the result is the same as in Song and Xu [15]. From the velocity profiles in Figure 4, we found that no matter whether CSF was included or not, the structure of time-averaged current profiles in the four cases were resembled. When including CSF, the magnitude of the surface current decreased, and the direction of the surface current rotated clockwise. Here, for a zonal wind, the surface stokes velocity can be approximated as U s0 = 0.24 m/s with the Stokes depth scale h s as 5 m [37,39]. The magnitude and direction of the surface current in the four cases are shown in Figure 4c, where in the cases of δ = 0 and δ = 0.6, surface currents rotate about 30°r ight to the wind and the velocity spirals are very flat, which is obviously different from the classic Ekman spiral, that is mainly due to the depth dependence of eddy viscosity in KPP. It is likely a convincing factor in explaining the discrepancies between observed surface currents and the classic Ekman theory. Once CSF is included in the Ekman model, this is equivalent to adding extra wave-induced stress to partly balance the wind-driven force. Therefore, relative to the cases without CSF (δ = 0 and δ = 0.6), a decrease in surface current was found and the surface flow rotated anticyclonically in the cases of δ = 0 + CSF and δ = 0.6 + CSF.   When including CSF in the Ekman model, the time-averaged surface current solid black line in right panel rotates clockwise relative to without CSF (dotted line in the left panel, and the magnitude of the surface current reduced, as part of the turbulent stress is used to balance the wave stress caused by CSF. It is interesting that the velocity vector traces a closed loop circle like the arabic numeral '8' when including CSF in the Ekman model (Figure 5b). This structure is related to the strength of CSF, which is directly proportional to the latitude and Stokes drift. As shown in Figure 6, under the same τ a , U s and δ, with the increase of the latitude, the CSF enhanced, which leads to a larger disturbance in a nonstationary Ekman model, and results in a changeable loop circle of surface current. This also explains the time-averaged surface current direction in the case of δ = 0.6 + CSF, which did not rotate anticyclonically to the case of δ = 0.6, but a little cyclonically, shown in Figure 5. Further insight into the Fourier series solution can be found in Figure 7, under the same forcing as shown in Figure 5. The tendency of the velocity varying along with time and depth is similar to the low-level jet stream in the atmosphere [6,21]. In the first half of 24 h, A v gradually reduces with time, and the Ekman layer depth starts from the deepest end. As A v decreases towards its midday minimum, the Ekman layer begins to shoal, and it can be easily found from the lower zero zonal velocity line from −h b to −0.6h b in Figure 5a, and the zero meridional velocity line form −0.6h b to −0.3h b in Figure 5b. A surface-intensified diurnal jet develops, associated with a high shear near the surface layer. Below the high-shear region, weak anticyclonic oscillations with an upwardpropagating phase are developed. Instead of this pattern, which is often due to inertial waves with downward energy propagation in the near-surface, here it is aroused by the diurnal cycle in eddy viscosity, while the right panel, including CSF, gives an enhancement of the shear and oscillation in the subsurface, along with a decrease in the magnitude of the surface current.
The zonal momentum balance is shown in Figure 8; the top three do not consider CSF, while the bottom three include the effect of CSF. It is similar to the inertial oscillations observed in the atmospheric boundary layer. Near the surface, there is an acceleration and deceleration of the flow on the side of the diurnal jet maximum. Deep in the layer, there are upward-propagation signals in acceleration, as a signature of inertial oscillations. This pattern is also like the results of Wenegrat and McPhaden [6] under a constant vertical eddy viscosity with the same time-dependence. After considering the CSF in the Ekman model, the basic pattern of momentum balance is unchangeable, but with an intensified up-propagation and vertical shear for the three momentum balance terms.

Rectification
Understanding the effect of time variability of eddy viscosity is important to reflect the true nature of Ekman dynamics. The colored dots in Figures 5 and 6 already reveal a local rectification effect of diurnal eddy viscosity, such as in Figures 7 and 8. Similarly to Wenegrat and McPhaden [6], they define the diurnal average of a variable X(t) as For comparison with the steady solution (δ = 0), which can be used to evaluate the rectification effects of the diurnal cycle in eddy viscosity, here, ω/2π is a period of the diurnal cycle which equals to 24 h. Additionally, with a simple normalized measure of rectification asX the bar notation represents the steady-state solution, assuming no time varies in eddy viscosity. The rectification effect of diurnal eddy viscosity relies on the strength of δ and its own vertical structure. For KPP shown in (4), that is the latitude (represented by f cor /ω) and u * . Under a constant wind stress force, two parameters f cor and δ determined the temporal and vertical structure of diurnal A v . Compared to Wenegrat and McPhaden [6], we have adopted a different depth-dependent eddy viscosity (KPP) and solve the Ekman problem in a Fourier series expansion, where a different rectification effect would be expected. The angle of the time-averaged surface current θ sur f as the function of f cor /ω and δ is shown in Figure 9, and the alteration of θ sur f due to diurnal varying A v is almost linear to f cor /ω when δ is smaller than 0.2, while with the increase of δ, the linear relationship with f cor /ω intensifies and the modification of θ sur f enlarges significantly. This is very different to the non-depth-dependent eddy viscosity shown in Wenegrat and McPhaden ([6], Figure 11). The friction velocity u * also has an influence on the structure of A v ; however, in ideal cases here, we did not consider the situation of time-varying wind and waves, so the rectification effect shown here is only arising from temporal variability of A v . Both the normalized rectification values of surface velocity (û R andv R ) and shear (û zR andv zR ) are shown in Figure 10. The velocity rectification intensifies with the increase of δ, when δ < 0.4, and the rectification of surface velocity is smaller than 0.1. The distribution ofû R andv R holds a decreasing tendency near the equator, probably due to the increase of absolute velocity as f cor → 0. However, the rectification of the shear is also proportional to δ and almost unrelated with the latitude, with robust growth at high δ, and the maximum values ofû zR 1 andv zR 4 as δ → 0.9. In general, the velocity shear is more strongly rectified than velocity, and the magnitude of rectification values are dominated by δ and weakly related to the latitude. In the middle and high latitudes, the rectification values in meridional are generally larger than in the zonal, as shown in Figure 10, wherev zR is obviously larger thanû zR , andv R is a little greater thanû R . It can be interpreted from the ellipse circles in Figure 11, that they project the major and minor axes of ellipse circles into the zonal (x) and meridional (y) directions, and we can get the corrected major and minor axes, which reflect the scale of variation in zonal and meridional velocity. Since the projected zonal axis is generally longer than the projected meridional axis, the oscillation frequency in the meridional axis would be higher in a diurnal period (24 h). This gives a proper explanation for the intensified rectification values in meridional velocity and shear.
To interpret the mean turbulent mixing due to rectification of diurnal A v , an effective eddy viscosity was defined to make the momentum balance equation without a tendency term satisfied for the mean velocity [1,2,6] A which is expressed in complex form, and there are two components of A v E f f , and any nonzero imaginary (or nonzero meridional) component in A v E f f stands for evidence of a rectification effect. The eddy viscosity in the KPP scheme holds a Gaussian vertical structure, and the normalized effective eddy viscosityÂ vE f f also holds the same convex profile shape, but the extents of magnitude and depth vary somewhat due to the rectification effect ( Figure 12). In the near surface,Â vE f f gets a little smaller thanÂ (black dashed line, which stands for non-time-dependent eddy viscosity), and strengthens the diurnal mixing δ, and the stronger decline inÂ vE f f . Below the near surface,Â vE f f gradually gets close toÂ and sequentially exceedsÂ at a different depth, where for δ = 0.9, the critical depth is about z/h b = −0.2, while for the weaker strength of δ, the deeper the critical depth is. Angles of the effective eddy viscosity θ A v E f f are generally positive under diurnal eddy viscosity and less than 45°, which means the effective viscosity rotates anticyclonically relative to the local mean shear. The angles change slowly with the depth, except near the surface and bottom layer. Strengthened diurnal mixing of δ means larger angles. Fluctuation of angles occurring at the bottom of the boundary layer is due to A vE f f → 0. The bottom two panels of Figure 12 give the magnitude of the real (zonal) and imaginary (meridional) components ofÂ vE f f , respectively. The intensified imaginary (meridional) component of A y vE f f at different δ interprets the growth of θ A v E f f in Figure 12b, which again confirms the existence of the rectification effect. It can be inferred that large δ deepens the penetration depth and amplifies mixing efficiency, and this generally follows from a previous study [1,2,6]. Latitudes also have an influence on the magnitude and angle of effective eddy viscosity, but do not revise the essential characteristics of A v E f f , since the rectification values of velocity and shear demonstrate a weak correlation with latitude.
It should be mentioned that the rectification effect shown above solely arises from the diurnal variation of eddy viscosity. The wave-modified terms in momentum Equation (1) are not considered in this section, for the evolution of wave properties will interfere with the rectification from time-varying eddy viscosity, and this issue will not be pursued further here.

Summary
Under the assumption of time-and depth-dependent eddy viscosity, we presented a Fourier series solution for a wave-modified Ekman model, and illustrated how the diurnal eddy viscosity alters the ocean's response to steady surface wind stress and rectifies time-averaged velocity.
In nature, the eddy viscosity does depend on time, although the dependence might be more complicated than we assumed. The time-related term in (2) breaks the structural constraint of eddy viscosity on time, and provides a simple insight into the physical realism of turbulent mixing. It has been shown that time-varying eddy viscosity modifies both the magnitude and vertical structure of ocean currents. Though diurnal period variation in eddy viscosity introduces no significant modification in time-mean velocity (Figure 4), an obvious upward propagation inertial oscillation appears in the vertical structure of velocity and momentum balance terms (Figures 7 and 8), along with a shoal Ekman boundary layer. These tendencies rectify the surface current in Ekman layer as a lowfrequency flow, like a low-level jet stream in the atmosphere [21].
Temporal dependence other than the vertical structure of eddy viscosity affects the rectification effect on velocity shear; therefore,û zR andv zR are highly linearly related to δ, but less related to latitude. Velocity shear is more strongly rectified than velocity, and rectification values in the meridional axis are larger than the zonal, probably due to the zonal steady wind force. Some of these results are similar to previous studies and some are new, particularly those related to rectification from time-and depth-dependent eddy viscosity in the KPP scheme. In this work, we did not attempt a comparison with observational data, and this study just aimed to emphasise a general lack of clarity regarding the scale of spatial and temporal averaging implicit in Ekman-type models.
Only the CSF was simply considered in the discussion of the solutions, as pointed out by Shrira and Almelah [7], where even under a constant wind, the Stokes drift cannot remain constant, and the dominant wavelength increases with time on a time-scale which might be comparable to or smaller than the scale of an Ekman layer formation. At present, there are very few studies of Ekman currents subjected to variable winds interacting with an evolving wave field. Rarely, wave-current-coupled numerical modeling can clearly illustrate the wave-current interaction. Since our Fourier series expansion solution can be extended to include wave effects, further work needs to be pursued for this issue.