Ekman layers in the Southern Ocean : spectral models and observations , vertical viscosity and boundary layer depth

Ekman layers in the Southern Ocean: spectral models and observations, vertical viscosity and boundary layer depth S. Elipot and S. T. Gille Scripps Institution of Oceanography, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0230, USA Scripps Institution of Oceanography and Department of Mechanical and Aerospace Engineering, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0230, USA now at: Proudman Oceanographic Laboratory, 6 Brownlow Street, Liverpool L3 5DA, UK Received: 8 December 2008 – Accepted: 19 December 2008 – Published: 12 February 2009 Correspondence to: S. Elipot (ship@pol.ac.uk) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
The Southern Ocean is believed to be a primary location of surface ocean mixing as a result of wind energy input, and this is of relevance for the global oceanic circulation (Wunsch and Ferrari, 2004).Large et al. (1997) stressed that observations of mixing processes from this region are needed to constrain general circulation models.A number of recent studies have evaluated mixing processes in the Southern Ocean, both in the deep ocean (e.g.Naveira Garabato et al., 2004;Sloyan, 2005) and in the upper ocean (e.g.Cisewski et al., 2005;Thompson et al., 2007).Nonetheless, we still lack observations of near-surface mixing on large scales.This study focuses on mixing processes that occur in the oceanic boundary layer (OBL) and that are linked to the wind-forced input of momentum to the upper ocean.
Most of our understanding of the ocean's response to wind forcing at the local scale has been framed in terms of Ekman (1905) theory, usually used to assess ocean response S. Elipot and S. T. Gille: Ekman layer in the Southern Ocean to constant or steady forcing.Steady conditions are rarely achieved in the real ocean, and, as a consequence, the steady Ekman spiral has proved difficult to observe.Only through extensive spatial and temporal averaging was it demonstrated to exist to some degree (e.g.Price et al., 1987;Wijffels et al., 1994;Chereskin, 1995).While predicted Ekman transports agree well with observations (e.g.Price et al., 1987;Chereskin and Roemmich, 1991;Chereskin, 1995;Wijffels et al., 1994;Schudlich and Price, 1998), predictions for the detailed vertical structure of the wind-driven velocities have been less satisfactory.Generally, an observed Ekman spiral appears more "flat" than the theoretical one derived from the "classic" steady model with a constant vertical viscosity K and an infinite ocean.This mismatch is an indication either that the velocity magnitude decays with depth more rapidly than the velocity vector rotates away from the wind stress direction or that the shear is predominantly downwind (Chereskin, 1995;Schudlich and Price, 1998;Price and Sundermeyer, 1999).K, which represents the "mixing", can be estimated by fitting observations either to the decay of speed with depth or to the velocity rotation at depth.Estimates obtained in these two ways can differ by an order of magnitude (Weller, 1981;Price et al., 1987;Chereskin, 1995;Lenn, 2006).Thus most studies have concluded that "Ekman theory" is unable to reproduce the observed detailed vertical structure of winddriven currents.However, Ekman's (1905) original theory included solutions for transient winds, and more recent work (e.g.Lewis and Belcher, 2004) has considered transient solutions with a range of different forms of vertical viscosity and upper ocean boundary conditions.
In this paper, rather than simply assuming that because Ekman theory does not perform well in the steady-state case it cannot apply under any circumstances, we systematically evaluate the behavior in the frequency domain of Ekmantype models.At sub-inertial frequencies, ocean observations indicate that turbulence closure models are good predictors of wind-driven velocities, while at the inertial frequency, slab-like models, which effectively have infinite vertical viscosity, are more successful (Weller and Plueddemann, 1996;Plueddemann and Farrar, 2006).We confirm here that the observed and theoretical ocean response to wind forcing strongly depends on time scale.Here we consider three different formulations for viscosity and three different formulations for the boundary condition at the base of the OBL.Together these produce nine different Ekman-type models, some of which have been investigated previously (e.g.Ekman, 1905;Gonella, 1972;Thomas, 1975;Madsen, 1977;Jordan and Baker, 1980;Lewis and Belcher, 2004), though we have found no previous comprehensive study of their frequency characteristics.The nine models each exhibit different behaviors in the frequency domain (see Appendix A).
Here we first identify the viscosity and boundary condition formulations that are best able to capture the observed relationship between time-varying wind and ocean velocity.Second, we ask whether the best of these time-varying Ekman-type models are sufficient to capture the full physics of wind energy input to the OBL.
The scrutiny that this paper gives to Ekman-type models might seem misplaced, since Ekman theory is not used in modern ocean general circulation models (OGCMs).Instead, modern OGCMs now typically use a turbulence closure boundary layer model such as KPP (Large et al., 1994), which includes a non-linear form for the vertical viscosity as well as a dependence on surface buoyancy fluxes and therefore must be solved numerically.In many regions, buoyancy forcing appears to play a significant role in determining the velocity structure of the upper ocean.For example, Price and Sundermeyer (1999) showed that the deepening and shoaling of the surface mixed layer by diurnal solar forcing could result in the time-mean spiral structure of wind-driven currents.Moreover, to model properly the wind-driven near-surface currents, local stratification should be taken into account, since it determines how wind-induced momentum penetrates through the water column (e.g.Plueddemann and Farrar, 2006).Ekman-type models do not explicitly represent buoyancy forcing or stratification, but that does not mean that the models are useless.We explore the possibility that the relevant effects of buoyancy forcing and stratification can be captured through an optimal choice of viscosity and boundary layer parameters.Moreover, Ekman models have other virtues.They are textbook classics with a long legacy, and they continue to inform our physical intuition about the upper ocean.They are analytically tractable, meaning that the influence of specific parameters can easily be explored.And simpler forms of Ekman-type models have been used in a number of recent studies of drifter data (e.g.Niiler and Paduan, 1995;Ralph and Niiler, 1999;Rio and Hernandez, 2003).
Looking at the characteristics of OBL models at different time scales comes down to considering their spectral characteristics as they appear in the transfer function with wind stress as input and ocean velocity as output (e.g.Gonella, 1972;Weller, 1981;Rudnick and Weller, 1993).Here, theoretical transfer functions are compared to the transfer functions estimated from surface drifter data from the Southern Ocean.The observed transfer functions are derived by carrying out cross-spectral analysis for surface drifter trajectories and wind stress interpolated onto drifter positions.
This paper is organized as follows: in Sect.2, the concept of a transfer function for vector input and output variables is interpreted in the context of OBL dynamics.In Sect.3, the mathematical steps leading to the transfer function expressions from the horizontal momentum balance equation are given.(The general characteristics of these transfer functions and their limiting behaviors are discussed in Appendix A. These functions can be graphically represented as a function of frequency and depth.)The oceanic and atmospheric datasets used in this study are described in Sect.4, and the methodology used to estimate the observed transfer functions in the Southern Ocean is given in Sect. 5.The results of fitting the modeled transfer functions to the observed ones are given in Sect.6 and a discussion of models' performance is found in Sect.7. Finally Sect.8 provides a summary.

Fourier series decomposition for a vector time series
A vector time series (here of the wind stress, drifter velocity or ageostrophic velocity) can be represented as a single complex Fourier series if it is treated as periodic with period T : where u and v are the zonal and meridional velocities, respectively; t is the time, and i= with the complex Fourier coefficient C k : Each component is a vector rotating with time.The hodographs for these vectors are counterclockwise-rotating circles for positive frequencies and clockwise-rotating circles for negative frequencies.For each rotary component, the absolute value of C k indicates its magnitude.

Theory of the transfer function for vectors
For our analysis, the local wind stress vector at the air-sea interface, τ (t), is interpreted as the input of a causal linear system.The output of this system is the ocean velocity vector u(t, z) at depth z.The velocity u(t, z) at time t can therefore be thought of as a convolution of the wind stress with the impulse response function h(t , z), where t is time lag, and z is depth (e.g.Bendat and Piersol, 1986, p. 189): Taking the Fourier transform 4), the convolution theorem linearizes the relationship: where U, H, and T are the Fourier transforms of u, h, and τ , respectively.At any given frequency ν, the transfer function H is complex valued.What is the interpretation of H? Assume the wind stress forcing is monochromatic (i.e. its Fourier series has only one non-zero component) with frequency ν 0 >0 and a magnitude of 1 N m −2 .Thus: The hodograph of such a wind stress is a counterclockwiserotating circle.Its Fourier transform can be defined with the help of the delta function, i.e.T(ν)=1×δ(ν −ν 0 ) (in units of N m −2 s).The resulting ocean velocity u(t, z) is the inverse Fourier transform of U(ν, z): = H(ν 0 , z) exp(+i2π ν 0 t).
(10) Thus, in this example, u(t, z) is a vector rotating with the wind stress at frequency ν 0 , and its Fourier series has only one non-zero component.The velocity vector has a constant deflection angle with respect to the stress vector, which is given by the phase of the complex number H(ν 0 , z) (in units of kg −1 m 2 s).If the rotating wind stress has a magnitude of 1 N m −2 (roughly equivalent to a 10-m wind speed of 20 m s −1 ; e.g.Large and Pond, 1981), then the absolute value of H(ν 0 , z) indicates the speed of the upper ocean currents.
In Appendix A, the theoretical and observed transfer functions are plotted in the complex plane.The axes can be thought as being fixed in a reference frame rotating with the wind stress vector, with the x-axis aligned with the wind stress vector.This representation is independent of the coordinate system, and it is particularly appropriate for studying the angular relationship between the wind-driven ocean velocity and the wind stress on global scales.This type of analysis is reminiscent of the averaging method developed by Price et al. (1987), where the signal-to-noise ratio of the wind-driven velocities is improved by projecting them into time-averaged along-and cross-wind directions.

Equation of motions
Our objective in this section is to evaluate how analytic Ekman-type models represent the relationship between wind and upper-ocean velocity.For a horizontally homogenous OBL, in the absence of pressure gradients the linearized horizontal momentum balance is: where u(t, z) is the horizontal velocity forced by the wind stress τ (t, 0), f the Coriolis parameter, and ρ the density of seawater.For consistency, the vertical coordinate z is taken positive downwards, and z=0 is the mean ocean-atmosphere interface.The "mixing" is written as a vertical flux of momentum per unit mass u w , where w is the vertical component of the velocity (positive downward).• represent the "fast" time average and primes the turbulent fluctuations that are typically not resolved by largescale oceanic observations.This flux defines a turbulent or Reynolds stress (per unit mass) acting on the large-scale circulation (e.g.Pedlosky, 1979): Following the concept that turbulent momentum fluxes are down-gradient and that they follow a Fickian law akin to what occurs at the molecular level, this turbulent stress is written as a turbulent coefficient K, the vertical viscosity, multiplied by the vertical shear of horizontal velocity: This parameterization provides a first order turbulence closure scheme of the Reynolds equations for the velocity in the OBL.It yields a linearized equation of motion conveniently written in terms of u only.Using Eq. ( 13), the momentum equation becomes: where K depends on depth only.The Ekman layer physics is governed by the vertical form of K and by the depth of the OBL.To obtain H, for each OBL model Eq. ( 14) is Fourier transformed to obtain an ordinary differential equation in z for U(ν, z): Then, using the Fourier transformed boundary conditions, a solution for U(ν, z) is found in the form given by Eq. ( 5) (see Sect. 3.3).
A number of studies have solved Eq. ( 14) explicitly for u(t, z) using a variety of vertical profiles for K(z) and applying several types of boundary conditions.For example, Lewis and Belcher's (2004) derivations of the timedependent solutions showed that if a constant wind-stress boundary condition is employed, then the lower boundary condition controls the damping scale, viscous or inertial, of the transient terms (in the form of inertial oscillations).Here, Eq. ( 14) is solved in the spectral domain.In many of the cases, our spectral solutions are modified versions of the time-mean terms of the solutions presented by earlier authors (Ekman, 1905;Thomas, 1975;Madsen, 1977;Jordan and Baker, 1980;Lewis and Belcher, 2004).

Parameterization of the vertical viscosity
We consider nine models arising from three different vertical profiles for K(z), and three different bottom boundary conditions.These models are sketched in Fig. 1.The model number (1, 2 or 3) designates the vertical profile of K and the letters (a, b or c) indicate the bottom boundary condition.Models 1a, 1b, and 1c have a constant viscosity K=K 0 (first row of Fig. 1), as in Ekman (1905).
Models 2a, 2b, and 2c have a viscosity that increases linearly with depth and that vanishes at the surface: K(z)=K 1 z (second row of Fig. 1).This linear increase in K with depth is physically justified, because it assumes that winddriven turbulent eddies are larger further from the surface, and therefore that the turbulent viscosity is larger at depth (e.g.Prandtl, 1952).For small z, a linear profile implies that the velocity should approximate a logarithmic profile as for a wall-bounded shear flow (e.g.Kundu and Cohen, 2002, p. 528), analogous to a linear K profile used for the atmospheric boundary layer (Tennekes, 1973).A similar profile has been predicted for the oceanic boundary layer (e.g.Madsen, 1977;Jordan and Baker, 1980;Thomas, 1975;Craig et al., 1993).
In models 2a, 2b, and 2c, zero values of K at the surface lead to a singularity.Models 3a, 3b, and 3c are designed to avoid this by using a viscosity that is finite at the surface and that increases with depth: K(z)=K 0 +K 1 z (third row of Fig. 1).The linear part of the viscosity profile is again justified by the mixing length argument.The constant K 0 allows the top boundary condition to be satisfied exactly without requiring approximations of the general solutions close to the surface.The addition of K 0 is equivalent to introducing a roughness length z 0 , so that the surface viscosity is K 0 =z 0 ×K 1 .Note that such a vertical profile for K approximates near the surface the cubic vertical profile implemented in the KPP model of Large et al. (1994).
Table 1.Mathematical expressions for the transfer functions H(ν, z).δ 1 = 2K 0 /(2πν+f ).δ 2 =K 1 /(2πν+f ).I n and K n are the nthorder modified Bessel functions of the first and second kind, respectively.For conciseness in the following table ζ (x)=2 i(z 0 +x)/δ 2 .The numbers on the left column designate the family of model and the letters on the first row designate the type of bottom boundary condition (see text).

K(z)
a -infinite layer b -one layer c -one and a half layer

Boundary conditions
For all models, the surface boundary condition matches surface wind stress to turbulent stress in the upper ocean.The boundary condition in the time domain and its corresponding Fourier transform are: This condition cannot be satisfied exactly when K vanishes at z=0 in models 2a, 2b, and 2c.Instead it is taken as a limit.
For the bottom boundary condition, three cases are considered: 1. Models 1a, 2a, and 3a are for a homogeneous ocean of infinite depth, and the corresponding bottom boundary condition specifies that the wind-driven velocity tends to zero: 2. Models 1b, 2b, and 3b are 1-layer models, with a homogeneous wind-driven finite layer of thickness h, at the bottom of which the velocity goes to zero: 3. Models 1c, 2c, and 3c are 1 and 1/2-layer models, consisting of a homogeneous wind-driven layer of thickness h, at the bottom of which the stress and hence the velocity shear go to zero, but non-zero velocity is still possible: (Price and Sundermeyer, 1999 used this bottom boundary condition to study the influence of stratification on Ekman layers.)

Mathematical expressions and graphical representations
The derivations of the transfer functions for models 1a, 1b, 1c, 2a, 2b, and 2c are omitted here because similar derivations have been published previously (e.g.Gonella, 1972;Thomas, 1975;Madsen, 1977;Weller, 1981;Lewis and Belcher, 2004).The transfer functions for models 3a, 3b, and 3c, to the best of our knowledge, are new results but their derivation is trivial1 .The mathematical expressions for the transfer functions of the models considered in this study are given in Table 1.These show that the ocean's response depends nonlinearly on the frequency ν of the forcing, the depth z, the latitude through the Coriolis parameter f , the water density ρ, and the vertical viscosity K.As indicated in Table 1 the depth scales for the transfer functions (δ 1 for models 1a, 1b, and 1c and δ 2 for models 2a, 2b, 2c, 3a, 3b, and 3c) depend on viscosity and frequency.
Appendix A provides further detail about the structure of the transfer functions and especially how each combination of vertical viscosity profiles and bottom boundary conditions leads to a different frequency response of the model (see Table A1).One interesting characteristic of these functions is their limiting behavior when the non-dimensional depths z/δ n , (n=1, 2) tend to zero.This situation occurs close to the surface and also when the angular frequency of the forcing approaches the inertial angular frequency −f .
In the Southern Hemisphere, f <0, and the inertial frequency is −f/2π >0.For cyclonic (ν≤0) and sub-inertial anticyclonic frequencies (0≤ν<−f/2π) all of the models indicate that the velocity is to the left of the wind stress at the surface and spirals downward anticylonically, while S. Elipot and S. T. Gille: Ekman layer in the Southern Ocean for supra-inertial anticyclonic frequencies (ν>−f/2π), the velocity is to the right of the wind stress at the surface and spirals cyclonically.The zero-frequency ν=0, or time-mean, velocity at the surface is consequently to the left of the mean stress direction.

Data
The Surface Velocity Program (SVP) (Siedler et al., 2001) and the ongoing Global Drifter Program (GDP) both provide horizontal velocity data from surface drifting buoys (drifters) on a global scale.A standard SVP drifter has a Holey-Sock drogue centered at 15-m depth, linked by a tether to a subsurface float and a surface float that radio-transmits its positions to the ARGOS satellite array at an uneven time rate, depending on satellite coverage and the drifter's setup (Sybrandy and Niiler, 1991;Niiler et al., 1995;Lumpkin and Pazos, 2007).The NOAA Atlantic Oceanographic and Meteorological Laboratory (AOML) Drifter Assembly Center processes the raw position data and interpolates them using a kriging procedure (Hansen and Poulain, 1996), resulting in a time series of position x(t) and velocity u d (t, x(t)) at six-hour intervals.
In principle, the drifter motions represent the currents averaged over the 6.1 m length of the drogue.Vertical shear of velocity has been observed over this lengthscale from vector measuring current meters mounted at the top and the bottom of the drogue (Niiler et al., 1995).Here shear information was not collected, and in our analysis we interpret the drifter velocities to be at the nominal 15 m depth.
In the Southern Ocean between 30 • S and 60 • S, 2839 independent SVP drogued drifter trajectories are available from November 1989 to May 2003.Undrogued drifter data were discarded.We identified 666 trajectories from drogued drifters that were at least 40 days long from October 1992, at the beginning of the TOPEX/Poseidon altimeter mission, to August 2002, the date when the ECMWF ERA-40 reanalysis ends (see below).Coastal areas are avoided by discarding the points of drifter trajectories for which a dynamic height relative to 3000 decibars from the 1 • gridded historical atlas data by Gouretski and Jancke (1998) could not be interpolated linearly.When divided in 40-day long segments that overlap by 20 days, these trajectories provide 10 387 time series segments, shown in Fig. 2.
These segments are further sorted in 2 • latitudinal bands according to their mean latitude (color-coded in Fig. 2).The number of segments per band is listed in Table 2.These numbers are used to evaluate the number of degrees of freedom for the spectral estimates, as explained in Appendix B2.
Figure 3a reveals the latitudinal biases, due to the decrease in data segments south of 44 • S.
In Fig. 3b, the longitudinal distribution of the data segments indicates that the drifters are primarily from the Atlantic and Indian sectors of the Southern Ocean.The temporal distribution of the data segments (Fig. 4) suggests that the observations are weighted more heavily toward the second half of the decade but show little seasonal bias.The drifter dataset is also further divided into an austral winter subdataset (5282 segments) and a summer subdataset (5105 segments) to study the seasonal variability.The austral winter is taken to correspond to the months of April through September and the austral summer to the months of October through March.The nominal month of a 40-day trajectory segment is chosen here as the month of the median date of the segment.
In order to obtain an estimate of the absolute geostrophic velocity component of the drifter velocities, two satellite altimetry datasets were combined.The anomalies u g were derived from "Archiving, Validation and Interpretation of Satellite Oceanographic" data that are produced by the Centre Localisation Satellite (AVISO).These provide high-resolution maps (1/3 • ×1/3 • Mercator grid) by merging TOPEX/Poseidon (T/P) and ERS-1 and -2 altimeter measurements, using an objective analysis method (Ducet et al., 2000).These maps are available at 7-day intervals implying a Nyquist frequency of 1/14 cycles per day (cpd) for the geostrophic part of the signal.We computed the velocity anomalies from the zonal and meridional gradients of the height anomalies.To these, a time-mean geostrophic velocity ūg was added, computed from the Gravity Recovery and Climate Experiment (GRACE) satellite-derived dynamic topography available on a global 1 • grid (Tapley et al., 2005).This mean geostrophic velocity field was interpolated 60 58 56 54 52 50 48 46 44 42 40 38 36 34 32  linearly in space, and the velocity anomaly maps were linearly interpolated in space and time, at all the drifter positions, to obtain the absolute geostrophic velocity u g + ūg at the surface every 6 h along the drifter trajectories.Time series of the ageostrophic velocity u at 15 m are then obtained as the drifter velocity minus the absolute geostrophic velocity at the surface: This neglects the geostrophic shear in the upper 15 m of the ocean.Expendable bathythermograph data in the Drake Passage indicate a geostrophic shear of less than 10 −3 s −1 in the upper 400 m (J.Sprintall, personal communication, 2006), yielding a potential maximum 1.5 cm s −1 geostrophic velocity difference between the surface and 15 m.This is of the same order as other sources of noise in this study.
For wind data, we use European Center for Medium-Range Weather Forecasts (ECMWF) ERA-40 Project reanalysis wind stresses (Simmons and Gibson, 2000) obtained from the Data Support Section of the Scientific Computing Division at the National Center for Atmospheric Research.The zonal and meridional wind stress components are available four times daily at the times 00:00, 06:00, 12:00 and 18:00 UTC.The values are instantaneous and are given as forecasts valid 6 h after the re-analysis time.The data are released on a Gaussian grid with resolution of 1.125 • longitude by roughly 1.125 • latitude.These gridded winds were linearly interpolated to the drifter positions to obtain contemporaneous six-hourly time series of wind stress τ (t).

Estimating the transfer function from the cross-spectra
We estimate the transfer function from observations using a spectral approach.The transfer functions discussed in Sect. 3 satisfy: where S τ u is the cross-spectral density function between the wind stress and the ocean velocity, and S τ τ is the autospectral density function of the wind stress.Here rotary power spectral density functions are estimated by the periodogram (e.g.Bendat and Piersol, 1986), for a finite number of frequency bands ν k : er of the legend, the linear dependences of the phase on frequency between 0 and 1 cpd converted to ime lags in hours (for the wind product with respect to the ocean velocity) are: 1. where • is the expected value operation over an ensemble of time series segments of length T and • * is the complex conjugate.X k is the finite Fourier transform of x: here computed using a standard Fast Fourier Transform algorithm.
The length of the time series segments considered here is T =40 days with a sampling interval t=0.25 day leading to N=160 points in time; thus the formal Nyquist frequency is 1/(2 t)=2 (cpd), although high frequencies are not wellresolved, and we focus our attention on frequencies well below the Nyquist frequency.The frequencies considered are ν k =k/T =k/(N t), positive for k=0, . . ., N/2 and negative for k=−N/2 + 1, . . ., −1.The frequency resolution is theoretically ν r =1/T = 0.025 cpd, but in reality it is 50% larger (0.0375 cpd), because we applied a Hanning window to reduce spectral side-lobe leakage (Harris, 1978).Since the data are ultimately sorted in 2 • latitudinal bands between 30 • S and 60 • S, this resolution is sufficient to resolve the smallest difference in the inertial frequency from one band to the next, except between the two southern-most bands.
The transfer function linking ocean velocities to wind stress is calculated from Eq. ( 21): using the data sorted in 2 • latitudinal bands.The zero frequency component is representative of the mean wind-driven currents at 15 m, and the phase at zero frequency is the mean angle over 40 days between the wind stress and either drifter or ageostrophic velocity at that depth.Table 2 lists χ (0) for 2 • latitudinal bands.At all latitudes χ (0) is greater for the drifter velocity than for the ageostrophic velocity most likely because of the oceanic eastward drift of the ACC flowing in the same direction as the atmospheric westerlies.The variation of the mean angle with latitude is one example of latitudinal variations in the transfer function (see below).

Correcting an unexplained constant time lag
Transfer functions of vector quantities are computed using rotary spectra (e.g.Mooers, 1973).Rotary spectra allow us to identify the angular separation between vector quantities but cannot distinguish differences in vector orientation from differences in temporal phasing.In each latitudinal band, we found that the phase of the transfer function depended linearly on frequency, which corresponds to a constant time lag between the wind stress and drifter data.First, in order to investigate if this lag was data-specific, several other types of wind products from the ERA-40 ECMWF Project re-analyses and the NCEP/NCAR Reanalysis Project (Kalnay et al., 1996) were tested.For the 52 • -54 • S latitudinal band, Fig. 5 shows the cross-spectral phases χ.Phases slope linearly with frequency for all products, but the slopes depend on the timing of the wind relative to the drifter measurements.This indicates that the time stamp of the data must be interpreted with care, particularly since wind products can be reported as instantaneous nowcasts, as forecasts (so that the time stamp precedes the actual wind by 6 h), or as time averages over 6-h intervals.In Fig. 5, the NCEP wind stress (black line), which is an average for the 6 h following the reported time, shows an expected constant time lag of 3 h with respect to the instantaneous ECMWF wind stress (red line), which is valid at the reported time.Surprisingly, the ECMWF winds show tilting phase lines (red line in Fig. 5) even when there is nominally no time separation between drifter and wind observations.A data-specific possible explanation for this is that there is a misalignment of the time stamps of wind and drifter data but why this lag should still be a function of latitude is unclear.We note also that Ekman models do not account explicitly for the presence of oceanic surface-gravity waves in the real ocean, and these could mediate the transfer of momentum from the atmosphere to the upper ocean.As such, a time lag could exist between the wind stress and the velocity at 15 m.How and why this lag would be a function of the frequency and latitude is unclear and further analysis, beyond the scope of this study, is required.
In conclusion, according to our analysis, it is found that the Ekman models considered here are unable to account for this time lag.Therefore, as a first step, we corrected the overall transfer function in each latitude band by removing the trend fitted in its phase between 0 and 1 cpd.Each trend correspond to a constant time lag reported in hours in Table 2.This lag decreases roughly from about 3 h to the north to about 1.5 h to the south of the region.

Influence of the wind slip
Surface drifters are excellent but not perfect water-followers, and their velocities contain an erroneous slip velocity caused by the direct action of the wind on the surface flotation buoy.Niiler et al. (1995) carried out experiments to measure wind slip in the tropical and Northeastern Pacific.They modeled the wind slip u s as: where w 10 is the 10-m wind velocity, R is the drag area ratio of the drogue to the other constituents of a drifter (40 for a SVP-type drifter), and a is a regression coefficient.Since no measurements in the field were obtained for winds stronger than 10 m s −1 , this model has not been validated for intense winds typical of the Southern Ocean: at drifter locations between 48 • S and 58 • S, the mean ECMWF reanalysis 10-m winds exceed 10 m s −1 , and the wind slip at these latitudes may be seriously underestimated (Niiler et al., 2003).
The standard wind slip in Eq. ( 26) was computed using ECMWF 10-m winds interpolated in time and space, and subtracted from the drifter velocities in order to obtain the wind slip-corrected velocities.Niiler et al. (1995) found that the best-fit values of a for either of two different types of drifters, TRISTAR or SVP Holey-Sock, were not statistically different.Their best estimate from the combined drifter datasets, a=4.63×10 −2 is used here.
We find that in general the wind slip correction reduces the magnitude of the real component of the transfer function, hence increasing the phase between stress and ocean velocity at all frequencies.The full consequences of this data modification are difficult to pin down, because the transfer functions and the optimization procedures are nonlinear.However, in general mean estimates of viscosity (see below) and boundary layer depth are not distinguishable within error bars from the estimates obtained when the wind slip correction is not applied.Furthermore, it does not make much sense to first remove a linear fraction of the wind in the form of an unvalidated wind slip correction and then subsequently to conduct a cross-spectral analysis between the "corrected" velocity and the wind stress.On the basis of these considerations, we have chosen here to present the results derived without the standard wind slip correction.

The cost function
Our aim is to determine the optimal boundary condition and vertical viscosity that best allow an Ekman-type model to represent the observed time-varying data.To do this, the observed transfer functions Ĥ are compared to the nine theoretical transfer functions H m listed in Table 1.We seek the optimal parameters for H m to minimize the cost function L, defined by the misfit between the observed and theoretical transfer functions: where | • | designates the absolute value.In the theoretical transfer functions, ρ is 1027 kg m −3 , the depth z is 15 m, and f is computed at the center of the 2 • latitudinal bands.
The L 1 -norm was selected rather than the L 2 -norm, because it performed better in the optimization procedure.Depending on the model considered, different algorithms were utilized for this nonlinear optimization.Details are given in Appendix B.
The weighting function, w(ν k ), is here the squared coherence γ 2 : and is estimated using Eq. ( 22).The normalized standard error of the cross-spectrum is theoretically inversely proportional to (γ 2 ) 1/2 (Bendat and Piersol, 1986), so that the best estimates of the cross-spectrum and hence of the transfer function are obtained when γ 2 is high and the weights w used to compute L penalize the frequency bands for which γ 2 is small.The minimum values of L resulting from the optimization procedures are plotted in Fig. 6d and dicussed in the next section.Near-surface data usually show that γ 2 is higher for anticyclonic frequencies than for cyclonic frequencies, and that it is higher at subinertial frequencies (Gonella, 1972;Weller, 1981;Daniault et al., 1985;McNally et al., 1989;Niiler and Paduan, 1995;Weller and Plueddemann, 1996;Rio and Hernandez, 2003;Elipot, 2006).Coherence is thought to decrease at lower and higher (absolute) frequencies mostly because of noise arising from other oceanic processes such as mesoscale geostrophic eddies or free inertial waves (Weller, 1981;McNally et al., 1989;Niiler and Paduan, 1995;Elipot, 2006).As a result of our choice of the weighting function w, our analysis of the data is most representative of sub-inertial motions and this should be kept in mind when interpreting the results.More specifically, the impact on the results of the zero frequency band or the frequency bands directly adjacent to this one is negligeable.While γ 2 will be reduced by noise, we find that we are able to reproduce the observed coherence fairly well a posteriori by using the theoretical expressions for H with the parameters estimated from the fitting procedure.Indeed, if one knows the transfer function H, the coherence can be predicted from the auto-spectra of the wind stress and the ocean velocity (e.g.Bendat and Piersol, 1986): While the transfer function H peaks at the inertial frequency (see Sect. 2), and the near-surface oceanic spectrum from drifter data has an approximate constant slope at subinertial frequencies (see e.g.Rio and Hernandez, 2003;Elipot, 2006), the wind stress spectrum slopes steeply at high frequencies (Gille, 2005;Elipot, 2006).Together these effects produce subinertial peaks for γ 2 , with higher coherence for anticyclonic frequencies (coherence squared level around 0.3) than for corresponding cyclonic frequencies (coherence squared level around 0.1).This readily translates into saying that up to 30% of the variance at sub-inertial anticyclonic frequencies and 10% at cyclonic frequencies are explained by the models.
6 Results of the fits

What are the best models for our observations?
To identify the optimal Ekman-model configuration, we assess which of the models has the smallest L, as plotted in Fig. 6d.We account for the uncertainty δL in this cost function, as defined in Appendix B4.Even with a quantitative cost function, no single model clearly outperforms all others at all latitudes.Figure 6a, b and c shows the viscosity coefficients K 0 and K 1 , and the boundary layer depth h, respectively, resulting from fitting the theoretical transfer functions of the models to the observed transfer functions in each 2 • latitudinal band.The error bars correspond to the mean absolute deviation from the mean of distributions drawn from a bootstrapping procedure (see Appendix B3).
Overall, the boundary condition c (no stress at the bottom or slip condition) is not helpful here.In all cases of vertical parameterization for K(z), the models with boundary condition c degenerate and are equivalent (see Fig. 6a, b and  Figs. A2, A1, A3) to the corresponding models with boundary condition a (infinite ocean): the optimal values for h are very large, ranging from physically acceptable for model 1c (O(10 3 m)) to unphysical and at the upper limit of the depth range explored by the optimization algorithms (see Fig. 6c).
Regardless of the representation of K(z), one-layer models (1b, 2b, and 3b) all perform significantly better than their counterparts with alternate boundary layers.
In summary, disregarding the "failing" models 1c, 2c, and 3c, the model performances are from best to worst: models 1b, 3b, 2b, 3a, 2a, and 1a.Revealingly, model 1a, the traditional Ekman model that has been tested extensively in previous studies, is the worst of these models.In the discussion  that follows, we focus first on the best model 1b.However, since this one returns parameters with unclear relationships to external environmental parameters (see Sect. 7), we also examine in detail the second best model 3b.

One-layer model with constant viscosity
Model 1b, with constant viscosity, a finite-depth boundary layer and a no-slip condition, should provide insight into the Ekman layer in the Southern Ocean.
Figure 8 shows the optimal parameters for this model for year-round data, as well as for summer and winter data.All 500 bootstrap estimates of each parameter are displayed in this figure.(See Appendix B1 for a discussion of the bootstrap procedure.)In some cases, the joint probability density function of K 0 and h (not shown) is bimodal rather than uni-modal, meaning that there are two distinct clusters of points in Fig. 8.This suggests that the subsampled data capture different types of oceanic conditions, while the scatter of each mode is intrinsic to random oceanic variability and random sampling of the data.
Throughout the Southern Ocean, this model indicates values for K 0 between 400×10 −4 m −2 s −1 and 1180×10 −4 m −2 s −1 (right panel of Fig. 8) and values for h between 30 and 50 m.The largest values of both K 0 and h are found between 40 • S and 50 • S. The joint distribution of bootstrap estimates of K 0 and h indicates a linear relationship between these two parameters: larger viscosities correspond to larger boundary layer thicknesses.This is consistent with the idea that K 0 represents turbulence stirred by the wind at the ocean surface, and h results from the same wind stirring.Linear fits between K 0 and h show that in most cases the minimum boundary layer depth is 15 m in the limit K 0 →0 since the optimization algorithm tries to force the drifter observations to be within the boundary layer.For this model h is found to be within a few meters of δ 1 (0)= √ 2K 0 /f , the exponential decay scale at zero frequency, which is the "depth of wind-currents" (divided by π) defined by Ekman (1905).
When the data are sorted by seasons, the scattering of the distributions increases, and at many latitudes the probability density functions of the bootstrap estimates indicate several modes (Fig. 8).However, the cost function is larger for the summer data than for the winter data (not shown), which makes the summer results less reliable.Thus the seasonal variability captured by this model is unclear.
Numerous studies have compared observed oceanic velocities with theoretical predictions from constant vertical viscosity models (see Huang, 1979;Santiago-Mandujano and Firing, 1990).Oceanic conditions, datasets, assumptions and processing methods all differ in these studies compared with our own.Broadly speaking, our results are consistent with those of Rio and Hernandez (2003), who also used surface drifter data and ECMWF wind stresses and who followed Ralph and Niiler (1999) in assuming a constant vertical viscosity within the Ekman layer.Rio and Hernandez (2003) filtered their data to retain a sub-inertial spectral band, and our cost function emphasizes the same frequency bands, so the similarities in our results are not surprising.Our viscosity estimates are however slightly larger and in fact closer to in situ estimates of about 10 −1 m −2 s −1 found near the Polar Front in the mixed layer in periods of strong winds (Cisewski et al., 2005).

One-layer model with linear viscosity with surface finite value
Model 3b has a linearly increasing viscosity with a finite non-zero value at the surface, K(z)=K 0 +K 1 z, and a finite boundary layer with a no-slip condition.The results and their seasonal variations are shown in Fig. 9.
This model degenerates to model 1b south of 50 • S since it returns values for K 1 that are not distinguishable from zero and values for K 0 and h that are not distinguishable within error bars from the values returned by model 1b.South of 50 • S few data are available and the transfer function estimates are more noisy.However, when the data are sorted by seasons, some of the bootstrap estimates, especially in summer, appear to continue the trend seen to the north of 50 • S.
North of 50 • S, the estimates of K 0 average (240±12)×10 −4 m −2 s −1 for year-round data, and they vary little with latitude.In contrast, h varies greatly with latitude.For year-round data north of 50 • S, h ranges between about 1400 m and 6000 m.It is smaller in summer compared to winter, and the latitudinal dependence is more pronounced in summer.In summer, h changes order of magnitude from north to south, increasing roughly from 350 m at 31 • S to 1925 m at 49 • S. In winter, h varies between about 2000 m and 6500 m, without clear latitudinal dependence.The implications of such large and unphysical values of h are discussed in the next section.Estimates of K 1 (lower left panel of Fig. 9) to the north of 50 • S range between 0.3 and 0.9 cm s −1 for year-round data.Two trends are noted for K 1 .First, for year-round data, it increases by a factor of 2.5 from north to south.Second, it increases from summer to winter by a factor 1.5 to the south and by 5.5 to the north.As discussed in the next section, the parameter K 1 is actually a friction velocity scale related to the wind stress.
Two-dimensional scatter plots of K 0 and K 1 of bootstrap estimates for each latitudinal bands and seasons (not shown) reveal a linear dependency between these two parameters.The larger K 0 is, the smaller K 1 .This is discussed in the next section.On the other hand, no relationship was found between h and either K 0 or K 1 .This suggest that the parameter h in this model captures a different signal in the data than do the K 0 or K 1 parameters.

Discussion
We are now left with two plausible models for the Ekman layer in the Southern Ocean, with two different parameterizations of the vertical viscosity.How do the parameters fitted for models 1b and 3b vary with respect to other environmental factors and what are their physical significance?

The relationship with the wind stress
The wind stress is the only forcing for Ekman models.Thus one might expect K and h to resemble the wind stress.geographical variability.The most noticeable feature in these two scales is that the seasonal variability disappears south of 48 • S.This is also the case for the viscosity scale u 2 * /f (not shown).
While model 1b provides the best match to the observed transfer functions, its optimal parameters h and K 0 show little of the latitudinal and seasonal variability that appears in the wind stress.This suggests that model 1b does not account for wind variability that should be important in the Ekman layer.Despite providing no simple dynamical insights, the optimal h and K 0 are within scaling ranges found in numerical studies of a neutrally stratified turbulent Ekman layer by Coleman et al. (1990).Our data show the latitudinally averaged ratio of h to u * /f for model 1b to be 0.32 for all data, 0.27 in winter and 0.45 in summer, comparable to the range 0.25-0.4found in numerical simulations.Similarly for model 1b, we found the average ratio of K 0 to the viscosity scale u 2 * /f to be 0.05 for all data, 0.04 in winter and 0.05 in summer, comparable to Coleman et al.'s range 0.03-0.08.
For model 3b optimal K 1 's and ECMWF u * 's are plotted in Fig. 10b.The coefficient K 1 , which has the units of a velocity, appears related to the wind stress.For models with linear viscosity, the linear coefficient is usually written K 1 =κu * (Thomas, 1975;Madsen, 1977), where κ is the Von Karman constant.Madsen (1977) assumed κ=0.4,but in the ocean or the atmosphere it is thought to be variable (Tennekes, 1973).To the extent that model 3b successfully captures oceanic variability, it gives us an unprecendented comparison between K 1 and u * .From our data (Fig. 10b), K 1 /u * =0.52 for all data, 0.64 in winter, and 0.33 in summer.In both seasons, this ratio increases with latitude.

The influence of stratification
When a slab layer model is used to simulate upper-ocean wind-driven velocity (Pollard and Millard, 1970) or to estimate the wind energy input to the mixed layer (D'Asaro, 1985;Alford, 2001), it is assumed that the wind momentum input is deposited uniformly throughout the wind-driven layer as a body force and this implies that the vertical profile of the wind-induced Reynolds or turbulent stress is linear.In these cases, the depth of the wind-driven "mixed-layer" is prescribed or limited, perhaps by a pre-existing stratification.The error bars for MLD, u * and u * /f are the standard error of the mean.The error bars for δ 2 (0) are obtained by formally propagating the errors from K 0 and K 1 taken as the mean of the absolute differences between the bootstrap estimates and the overall most probable estimate.
system through a linear drag term that is intended to represent radiation of energy out of the mixed/wind-driven layer.The drag coefficient is typically tuned to match the velocity observations, but it has been shown that this typically overestimates the wind energy input (Plueddemann and Farrar, 2006).In contrast, boundary layer models, such as KPP, that explicitly incorporate buoyancy forcing deposit momentum to a "surface layer" or shallowest layer.Then, the depth over which the vertical viscosity is enhanced by the wind momentum input, the boundary-layer depth (BLD), is usually deeper and is diagnosed by a criterion based on a bulk Richardson number relative to the top most layer of the numerical model.The simple idea is that the stratification limits the vertical penetration of turbulent momentum.However, in the tropical Pacific, Zhang and Zebiak (2002) found that KPP produced more realistic velocities when it was modified to deposit wind momentum as a body force over the whole BLD.
In "Ekman" models the wind-induced stress is a non-linear function of depth, and it is not associated with a constant body force per unit mass.In that case, energy is removed from the system only by dissipation through the shear induced stress and the downward radiation of energy by internal waves or the deepening of the mixed-layer is not mod-eled.This is clearly a limitation when modeling the real ocean (e.g.Plueddemann and Farrar, 2006).Our optimization procedure requires only that the BLD be less than 10 4 m, and the optimal BLDs obtained for model 3b are O(10 3 m), values that can at times exceed the water depth and that are clearly unphysical.One possible explanation is that such models are unable to extract enough energy from the system, and they set the boundary layer to be extremely deep to accomodate large wind energy input.
The KPP formulation uses a cubic profile for K(z).In a coarse resolution OGCM for the Southern Ocean, Large et al. (1997) found that the monthly-mean mixed-layer depth (MLD) and BLDs determined by KPP were comparable.However, on much shorter time scales when the stirring by the wind is intense, they noted that the BLD could be much greater than a MLD defined as an isothermal layer.For this study we compared the observed transfer functions to transfer functions derived from a KPP-style cubic profile of the vertical viscosity.This required numerical solution.The resulting viscosity estimates were indistinguishable from the estimates obtained by the linear viscosity models, because our estimated BLD was again unphysically large O(10 4 m), and the cubic profile approximated a linear profile near the surface, much like models 3a, 3b, and 3c.vertical viscosity did not produce plausible BLDs, it is clear that buoyancy fluxes play a role in the ocean boundary layer that Ekman models are unable to represent.
To investigate this further, we explored whether we could detect the influence of the stratification in our results.The climatological MLD determined from density profiles (Dong et al., 2008) was interpolated in space and time to the drifter positions.Mean values are plotted in Fig. 10a, as a function of latitude and season.MLDs and BLD h from the drifter data differ by an order of magnitude, with MLD being O(100 m) and h O(1000 m).Nonetheless, both exhibit common latitudinal and seasonal trends, implying that the stratification represented by MLD can be associated with BLD.
Interestingly, at each latitude band, the depth scale δ 2 at zero frequency (filled symbols in Fig. 10a) is close to the mean value of the MLD.This correspondence is found not only for year-round data but also for seasonally sorted data.
Whereas δ 1 (0) for models 1a, 1b, and 1c is a familiar scale of exponential decay, δ 2 (0) appears in a complicated manner in the expression of the transfer function for model 3b (see Table 1).We computed the ratio of the absolute value of the transfer function at the depth z=δ 2 (0) to the surface value, which is also the ratio of the velocity magnitudes at the same depths, using the optimum parameters.At the depth z=δ 2 (0), the current speed is about 15% of its surface value at 50 • S.This percentage increases to about 32% at 31 • S.These results imply that the shear is large and velocities greatly reduced at the "Ekman depth".
Overall, the model 3b results suggest that the wind-driven velocities penetrate deeper than the depth of the mixed layer but that the mixed-layer depth nevertheless controls the Ekman scale of the model.

Speculation about the sea surface roughness
The atmospheric boundary layer and the oceanic boundary layer interact with each other and create roughness along their interface (e.g.Melville, 1977).For the ocean, the roughness length z 0 is expected to be representative of the thickness of an unresolved, wave-enhanced sub-layer (Craig and Banner, 1994), just below the surface.Possible scalings for z 0 found in the literature include some multiple of u 2 * /g (Charnock, 1955) where g is the gravitational acceleration, the wavelength of the waves (Craig and Banner, 1994), or the significant wave height (e.g.Terray et al., 1996).The length z 0 needs to be considered in order to model correctly the vertical velocity profile as one approaches the boundary.For models 3a, 3b, and 3c, the optimization procedure was set up to conduct a search of the two parameters K 0 and K 1 , which were assumed to be independent.A scatter plot (not shown) of all bootstrap estimates of K 0 versus K 1 in each latitudinal band shows that they are actually linearly dependent; Fig. 11 shows the linear coefficient or roughness length z 0 =K 0 /K 1 for model 3b.
The roughness parameter is larger in the austral summer than it is in the austral winter, which is mostly a consequence of the seasonal variations of K 1 .An examination of Fig. 11 suggests no clear relationship between z 0 and MLD, wind stress, or the Coriolis parameter.Further investigation is required to link these estimates to oceanic conditions.

Summary
This paper has studied the frequency response of the ocean boundary layer to wind stress forcing.We used a series of Ekman-type models, so named because no explicit buoyancy forcing is considered and the turbulent vertical flux of horizontal momentum is parameterized by a first-order turbulence closure as first introduced by Ekman (1905) for the ocean.Models of this type are highly idealized in some respects, which might make them seem inappropriate for real ocean applications, but they have a long tradition in physical oceanographic literature and they are mathematically tractable, making them a natural starting point for any consideration of upper ocean physics.
Ocean Sci., 5, 115-139, 2009 www.ocean-sci.net/5/115/2009/ S. Elipot and S. T. Gille: Ekman layer in the Southern Ocean 131 We have sought a formulation of Ekman-type models to best represent the frequency response of upper ocean currents to time-varying winds.Three vertical profiles for the vertical viscosity are considered: a constant profile, a linear profile increasing with depth from zero at the surface and a linear profile increasing with depth from a finite value at the surface.Three boundary conditions for the bottom of the oceanic boundary layer are considered: an infinite depth layer with vanishing velocity at infinite depth, a finite depth layer at the bottom of which the velocity vanishes and a finite depth layer at the bottom of which the stress vanishes.Together these imply nine different models.The frequency response of each of these model is described by a depthdependent transfer function H.At each frequency, the phase of H gives the deflection angle of the oceanic velocity with respect to the instantaneous wind stress, and the magnitude of H indicates the magnitude of the oceanic velocity for a 1 N m −2 wind stress.
Parameters for the theoretical transfer function H are tuned to find the best match to transfer functions derived from Southern Ocean drifter observations, altimetry, and wind fields, and the success of the models has been evaluated.Results show that the classical Ekman model, with constant vertical viscosity and infinite depth, is among the least successful representations of the OBL.The model can be improved by treating the upper ocean as a 1-layer system and/or by allowing the vertical viscosity to vary with depth.
The best model to describe the frequency response of Southern Ocean drifter velocities to wind stress forcing is a one-layer model with a constant vertical viscosity.From 60 • S to 50 • S, the boundary layer is shallow, of O(30-35) m, and the viscosity averages 724×10 −4 m −2 s −1 with small seasonal variations of the order of ±15%.These latitudes correspond to the largest zonally-averaged wind stress values in the Southern Ocean with little seasonal variations.From 50 • S to 40 • S, the boundary layer is best described by a slightly deeper layer O(45-50) m, with associated increased constant vertical viscosity reaching over 1000 m 2 s −1 , however with very little seasonal variability.From 40 • S to 30 • S, the boundary layer is shallower again, O(35) m, and the viscosity is smaller, averaging 474×10 −4 m −2 s −1 .
For the latitude range from 30 • S to 50 • S, an alternate description of the Ekman layer is given by a one-layer model with a vertical viscosity that increases linearly with depth from a finite value at the surface.In this model, the boundary layer parameter and the vertical viscosity coefficient K 1 both appear to vary with wind.The boundary layer is much deeper than the mixed layer, with deepest values in winter and at latitudes where the wind is strongest.The vertical viscosity coefficient K 1 is O(10 −3 −10 −2 m s −1 ) and scales like the friction velocity, showing similar seasonal and latitudinal variations.The viscosity at the surface K 0 ranges between 10 −2 and 4×10 −2 m 2 s −1 and does not show obvious dependence on latitude, wind stress or MLD.The boundary layer depth parameter is O(10 3 m) and can exceed the ocean depth, implying that it has limited physical meaning.In contrast, the time-mean Ekman depth scale, K 1 /f does appear physically meaningful: it is close to the seasonally and latitudinally varying climatological MLD.
The models presented here are not able to explain fully the observed locally wind-driven variability in the upper ocean.Buoyancy fluxes, which are omitted from Ekman-type models, need to be considered to account for the unexplained portions of the upper ocean response to wind.

Limiting behavior of the transfer functions
The frequency and depth dependence of the transfer functions can be illustrated graphically.Figure A2 shows the transfer functions for models 1a, 1b, 1c, Fig. A1 for models 2a, 2b, 2c, and Fig. A3 for models 3a, 3b, 3c.These transfer functions are evaluated with numerical values for the viscosity K and the boundary layer depth h, chosen as optimal parameter fits for Southern Ocean observations (see Sect. 6) in the 40-42 • S latitudinal band.The plots shown here are representative examples of the zonally-averaged OBL in the Southern Ocean.Frequencies are plotted from −2 cycles per day (cpd) to 2 cpd, since the 6-h data have a Nyquist frequency of 2 cpd.The vertical variation of the transfer function is plotted as a line, color-coded by frequency.Each curve in these figures is analogous to the velocity hodograph as a function of depth, or what could be called an Ekman "spiral".The colored dots (on the lines in Fig. A2 or projected on the (x, y, h) plane in Figs.A1 and A3) give the transfer functions at 15 m for each frequency band.The observed transfer function at 15 m estimated from the data in the 40-42 • S latitudinal band is plotted on the (x, y) plane in the lower-right panels of Figs.A2, A1, and A3.For models 1a, 1b, 1c and 3a, 3b, 3c the transfer function at the surface as a function of frequency are plotted with gray curves.For model 1c, the transfer function at the bottom of the boundary layer is also drawn (lower-left panel of Fig. A2).

A1 Constant eddy viscosity models
For K=K 0 (models 1a, 1b, and 1c), the general solution of Eq. ( 15) is The steady case for model 1a is obtained from the expression in Table 1 by setting ν=0.This gives the "classic" timeinvariant Ekman spiral solution: where is the exponential decay scale.D E =π|δ e | is the "Depth of Wind-currents" defined by Ekman (1905), which is the depth at which the velocity is opposite in direction to the velocity at the surface.At non-zero frequencies, the exponential decay scale is modified and we define a frequency-dependent "Ekman depth": ) |δ 1 | represents the penetration depth of the wind-driven currents, which increases with the square root of K 0 , since a larger viscosity is expected to be representative of more vigorous turbulence, and is inversely proportional to the square root of the "wind rotation" ν * =2πν+f (Crawford and Large, 1996).Frequency ν * is a measure of the relative rotation in the local reference frame at the cyclonic frequency f/2π (units of s −1 ).When the frequency is inertial (ν=−f/2π), |δ 1 | goes to infinity.
The transfer functions for models 1a, 1b, 1c (first row of Table 1) are written in a way that emphasizes the angular separation at the surface.Table 1 shows that model 1a has an angular separation at the surface of ±π/4 for all frequencies, and it increases with depth, anticyclonically for sub-inertial frequencies and cyclonically for supra-inertial frequencies.For models 1b and 1c, the deflection angle is influenced by the finite thickness h of the boundary layer and can therefore differ substantially from π/4 at the surface.
We examine the behavior of the transfer functions near the inertial frequency, in the limit where 2πν→−f .For model 1a, the velocity at all depths is predicted to be nearly oriented at ±π/4 from the wind stress (see Table A1), and the magnitude of the response has an unbounded resonance.Model 1b and model 1c near-inertial behaviors are very different (see Table A1), and this emphasizes that choosing the right bottom boundary condition is potentially crucial for modeling high-frequency wind-driven currents.For model 1b, the inertial resonance is finite and downwind at all depths, and the vertical shear is constant.The inertial surface drift scales like h and inversely with K 0 .In contrast, for model 1c, the inertial resonance is infinite, the shear is zero, and velocities at all depths are at right angles to the wind direction.The transfer function scales inversely to h and is independent of the viscosity.This is an inertial slab-like behavior but since the shear is zero, there is no dissipation term to remove energy from the system.This forced inertial "mode" of motion is unlikely to represent real oceanic processes.Similarly, Lewis and Belcher (2004) found in the time dependent solution for model 1c that an undamped mode oscillating at the inertial frequency is excited when an impulsive stress is imposed on an ocean originally at rest, and they consequently abandoned this model as being unphysical.In Sect.6, we find that this model performs poorly, most likely because the data indicate a downwind inertial response.

A2 Linear viscosity models
For K=K 1 z (models 2a, 2b, 2c), the general solution of Eq. ( 15) is: where I n and K n are the nth-order modified Bessel functions of the first and second kind, respectively, and ) is a new frequency-dependent Ekman depth for models 2a, 2b, and 2c (and also for models 3a, 3b, 3c) that goes to infinity at the inertial frequency.A(ν) and B(ν) are determined by the boundary conditions.The surface boundary condition Eq. ( 16) is taken as the limit using first-order approximations  for the derivatives of the Bessel functions (Madsen, 1977).The mathematical expressions of the transfer functions for these models are given in the second row of Table 1 for the three bottom boundary conditions.Madsen (1977) and Lewis and Belcher (2004) both derived the transfer function for model 2a in Laplace transform form and inverted it to obtain the time dependent solution in the oceanic boundary layer.
The behaviors as z/δ 2 →0 are summarized in Table A1.These are obtained by retaining the first term of a series expansion for K 0 around 0 (see Table A2).
For model 2a, the imaginary part of the transfer function (the crosswind component of velocity) tends to a constant, while the real part (the downwind velocity component) is logarithmic and eventually goes to infinity.Model 2b presents a rather different limiting behavior than model 2a: Table A2.Limiting behaviors for small argument of the zeroth and first orders modified Bessel functions of the first and second kinds.is the Euler constant.
it predicts that near the surface, the oceanic boundary layer behaves like a logarithmic layer and that there is no crosswind component for the inertial response.The limiting behavior of model 2c is a combination of the limiting behavior of model 2a and model 1c: it has a logarithmic downwind component with a constant cross-wind component and also includes an "inertial" mode at right angles to the wind that is independent of the viscosity but dependent on the boundary layer depth.In Sect.6 we find that that model 2c fails in the sense that fitted values for h are physically too large.
For models 2a, 2b, 2c the singularity at z=0 is inconvenient, because the surface velocity is not defined.In order to obtain this surface "drift", Madsen (1977) evaluated the velocity at a depth z 0 from the theoretical surface.This distance is called the roughness length and for the case of an OBL could correspond to an unresolved sub-layer just beneath the surface where turbulence caused by waves (breaking or not) occurs.The size of z 0 is subject to much debate (e.g.Stips et al., 2005).Reviewing field and laboratory experiments, Madsen (1977) used a length of O(10 −2 m) and found that only the order of magnitude was relevant since a multiplicative factor of 2 for z 0 changed the surface drift magnitude and angle by only 10%.In Sect.6, we find that the fitted values for the linear coefficient K 1 in the Southern Ocean are one to two orders of magnitude larger than those

Fig. 2 .Fig. 2 .
Fig. 2. Drifter trajectory segments used in this study between 30 • S and 60 • S. The 40-day segments are colored according to their mean latitude, following a repeated 5-class qualitative colormap to distinguish one 2 • latitudinal band from the next.

Fig. 3 .
Fig. 3. (a) Latitudinal distribution and (b) longitudinal distribution of the median dates of the of the 20-day overlapping 40-day drifter trajectory segments.
Fig. 4. (a) Month distribution of the median dates, and (b) year distribution of the mean latitude of th overlapping 40-day drifter trajectory segments used in this study.
S. Elipot and S. T. Gille: Ekman layer in the Southern Ocean τ −6h NCEP τ average of NECP and NCEP τ −6h NCEP 10−m wind ECMWF 10−m wind ECMWF τ −6h hase of the cross-spectrum between the drifter ageostrophic velocities and various wind and wind for the data in the 52 • -54 • S latitudinal band.ECMWF stress, ECMWF 10-m wind and NCEP 10-m instantaneous values valid at the drifter time.NCEP stress −6 h is the average value valid over the h before the drifter time.NCEP stress is the average value valid over the next 6 h starting from the e.Average NCEP stress is the arithmetic average of these last two values.ECMWF stress −6 h is taneous stress value valid 6 h before the drifter time.A positive phase means that the ocean velocity eft of the wind.A positive linear slope of the phase indicates that the wind lags the ocean velocities.

Fig. 7 .Fig. 7 .
Fig. 7. Average cost function values for all data (left panel), data south of 50 • S (middle panel) and data to the north of 50 • S (right panel).Error bars are plus or minus the latitudinally averaged standard error for the cost function as derived in Appendix B.
For a stable planetary boundary layer, the relevant planetary scale is u * /f , where u * = √ |τ |/ρ is the friction velocity scale.Figure 10a shows u * /f , and Fig. 10b shows u * derived from the ECMWF wind stress.Since these scales are evaluated from the mean of the values of wind stress interpolated at the drifter locations, they should reflect the same seasonal and h and vertical viscosity K 0 for model 1b.The results for year-r k, for summer in red and for winter in blue.The overall optimum parameters are mbols and the bootstrap distributions are plotted with colored dots.

Fig. 8 .
Fig. 8. Boundary layer depth h and vertical viscosity K 0 for model 1b.The results for year-round data are plotted in black, for summer in red and for winter in blue.The overall optimum parameters are plotted with white-filled symbols and the bootstrap distributions are plotted with colored dots.

Fig. 9 .
Fig. 9. Boundary layer depth h and vertical viscosity coefficients for model 3b.The results for year-round are plotted in black, for summer in red and for winter in blue.The overall optimum parameters are plotted with white-filled symbols and the bootstrap distributions are plotted with colored dots.

Fig. 9 .
Fig. 9. Boundary layer depth h and vertical viscosity coefficients for model 3b.The results for year-round are plotted in black, for summer in red and for winter in blue.The overall optimum parameters are plotted with white-filled symbols and the bootstrap distributions are plotted with colored dots.

Fig. 11 .
Fig.11.Sea surface roughness estimates z 0 =K 0 /K 1 in 2 • latitudinal bands for models 3b.Note that no overall optimum estimates cannot be obtained south of 50 • S since the optimum K 1 ≈0.
Fig. A1.Transfer functions for model 1a with K 0 =574×10 −4 m −2 s −1 ; model 1b withK 0 =106×10 −4 m −2 s −1 and h=51 m; model 1c with K 0 =558×10 −4 m −2 s −1 and h=1528 m; f =−0.95×10 −4 s −1 corresponding to 41 • Sand an inertial frequency of approximately 1.3 cpd.Each curve is the transfer function as a function of depth for frequencies ν=−1.95 . . .1.95 cpd at 0.05 cpd interval, with lines color-coded by frequency.The black curves are the transfer functions at the zero-frequency.The transfer function at 15 m is indicated by a colored dot on each curve for each model.The gray curve joins the z=0 m points for all frequencies for models 1a, 1b and 1c.For model 1c a gray curve also joins the z=h points.The dotted lines indicate the x and y axes and the ± 45 • directions.The lower-right panel is the observed transfer function at 15 m in the 41 • S zonal band.

Fig. A1 .
Fig. A1.Transfer functions for model 1a with K 0 =574×10 −4 m −2 s −1 ; model 1b with K 0 =106×10 −4 m −2 s −1 and h=51 m; model 1c with K 0 =558×10 −4 m −2 s −1 and h=1528 m; f =−0.95×10 −4 s −1 corresponding to 41 • S and an inertial frequency of approximately 1.3 cpd.Each curve is the transfer function as a function of depth for frequencies ν=−1.95 . . .1.95 cpd at 0.05 cpd interval, with lines color-coded by frequency.The black curves are the transfer functions at the zero-frequency.The transfer function at 15 m is indicated by a colored dot on each curve for each model.The gray curve joins the z=0 m points for all frequencies for models 1a, 1b and 1c.For model 1c a gray curve also joins the z=h points.The dotted lines indicate the x and y axes and the ±45 • directions.The lower-right panel is the observed transfer function at 15 m in the 41 • S zonal band.
Fig. A2.Transfer functions for model 2a with K 1 =0.77×10 −2 s −1 ; model 2b with K 1 =0.42×10 −2 s −1 and h=56 m; model 2c with K 0 =0.77×10 −2 s −1 and h O(10 4 ) m; f =−0.95×10 −4 s −1 corresponding to 41 • S and an inertial frequency of approximately 1.3 cpd.See also the caption for Fig.A1.The theoretical transfer functions at 15 m depth are projected on the plane coinciding with the bottom of the axes.The real part of the transfer functions at ν=0 is projected on the (x, z) plane and the imaginary part on the (y, z) plane and these curves are drawn in black.Since these transfer functions are not defined at the surface, the curves curves start at the depth z=−0.1 m.The lower-right panel is the observed transfer function at 15 m in the 41 • S zonal band.
Phase of the cross-spectrum between the drifter ageostrophic velocities and various wind and wind stress data for the data in the 52 • -54 • S latitudinal band.ECMWF stress, ECMWF 10-m wind and NCEP 10-m wind are instantaneous values valid at the drifter time.NCEP stress −6 h is the average value valid over the previous 6 h before the drifter time.NCEP stress is the average value valid over the next 6 h starting from the drifter time.Average NCEP stress is the arithmetic average of these last two values.ECMWF stress −6 h is the instantaneous stress value valid 6 h before the drifter time.A positive phase means that the ocean velocity is to the left of the wind.A positive linear slope of the phase indicates that the wind lags the ocean velocities.In the order of the legend, the linear dependences of the phase on frequency between 0 and 1 cpd converted to constant time lags in hours (for the wind product with respect to the ocean velocity) are: 1.62, 4.69, −1.36, 1.68, 1.77, 1.26, 7.64.

Table 2 .
Characteristics of trajectory segments per 2 • latitudinal band.The lag is discussed in Sect. 5. χ(0) is the mean angle between the wind stress and drifter velocities u d or ageostrophic velocities u.

Table A1 .
Mathematical expressions for the limiting behaviors of the transfer functions for z/δ n → 0. δ 1