Turbulent universality and the drift velocity at the interface between two homogeneous fluids

The drift velocity U0 at the interface between two homogeneous turbulent fluids of arbitrary relative densities in differential mean motion is considered. It is shown that an analytical expression for U0 follows from the classical scaling for these flows when the scaling is supplemented by standard turbulent universality and symmetry assumptions. This predicted U0 is the weighted mean of the free-stream velocities in each fluid, where the weighting factors are the square roots of the densities of the two fluids, normalized by their sum. For fluids of nearly equal densities, this weighted mean reduces to the simple mean of the free-stream velocities. For fluids of two widely differing densities, such as air overlying water, the result gives U0 ≈ αV∞, where α 1 is the square root of the ratio of the fluid densities, V∞ is the free-stream velocity of the overlying fluid, and the denser fluid is assumed nearly stationary. Comparisons with two classical laboratory experiments for fluids in these two limits and with previous numerical simulations of flow near a gas-liquid interface provide specific illustrations of the result. Solutions of a classical analytical model formulated to reproduce the air-water laboratory flow reveal compensating departures from the universality prediction, of order 15% in α, including a correction that is logarithmic in the ratio of dimensionless air and water roughness lengths. Solutions reproducing the numerical simulations illustrate that the logarithmic correction can arise from asymmetry in the dimensionless laminar viscous sublayers.


I. INTRODUCTION
The flow adjacent to an interface between two homogeneous, turbulent fluids is a classical topic in fluid dynamics, which continues to stimulate current research because of its importance in a wide range of environmental and engineering contexts [1][2][3][4][5][6][7] . One aspect of this problem is the determination of the drift velocity: the mean interface-parallel motion at the interface. It is shown here (Sec. II) that a modest extension of classical scaling for these flows, supplemented by standard turbulent universality and symmetry assumptions, yields a quantitative, analytical prediction of this drift velocity for the general case of two homogeneous fluids of arbitrary relative densities (Fig. 1). Despite the classical nature of the derivation, this result does not seem to have been previously noted.
The main purpose of this note is to record this result and make it available in the literature. A second purpose is to explore briefly the nominal accuracy of the prediction and some possible sources and anticipated magnitudes of departure from it. While the analytical result for the drift velocity is exact under the universality and symmetry assumptions, it is important to recognize that those assumptions will be violated to some degree in any physical fluid system, and consequently the result must be understood only to provide a starting point for further, more complete analysis. This is particularly true for flows supporting significant interfacial wave activity, the complex effects of which on flow near the interface are themselves a focus of much current research 2,[8][9][10][11][12][13][14][15] . Nonetheless, the universality prediction suggestively reproduces a long-standing empirical rule for ocean wind drift that was reported already more than half a century ago 16,17 and is still in operational use 18 , according to which the surface drift current speed is approximately 3% of the wind speed.
Consistent with the classical nature of the derivation, the result is compared with two classical laboratory experiments 19,20 that include explicit measurements of the velocity profiles on both sides of the interface and are only modestly affected by interfacial waves. For the air-water interface, the characteristic profile of air velocity (perhaps above a wave boundary layer 21,22 ) is generally regarded as well-enough known 23 that more recent laboratory studies generally measure the velocity profile only in the water, and primarily focus instead on wavestate effects 24,25 . The result is compared also with a pioneering numerical study of coupled gas-liquid flow, which is conducted in a symmetric, dimensionless setting that anticipates the development by which the drift velocity prediction is obtained 26 .

II. UNIVERSALITY RESULT
Consider two homogeneous fluid layers, the first with density ρ 1 overlying the second with density ρ 2 > ρ 1 , in mean differential motion, with turbulent boundary layers adjacent to the interface. Suppose that the flow in both layers is statistically homogeneous in the horizontal and stationary in time, and let the constant, horizontal free-stream velocities outside the upper and lower fluid boundary layers be denoted by V ∞ and U ∞ , respectively.
In general, the interface will be distorted by interfacial waves induced by the turbulence and the shear layer. Let z be the quasi-vertical distance from the interface in a curvilinear coordinate system, so that z = 0 defines the instantaneous position of the distorted interface.
Under the assumptions of homogeneity and stationarity, the Reynolds' or ensemble averaged velocities in the upper and lower layers can then depend only on z, and may be denoted respectively by V(z) and U(z). Note that the implicit use of an interface-relative coordinate implies restrictions on the degree of distortion and connectedness of the interface, which are discussed further in Section IV.
At the interface, standard conditions 27 of continuity of velocity and stress are required to hold for the instantaneous flow, and so must also hold for the averaged flow: where the stress τ (z) is the mean turbulent vertical flux of horizontal momentum and the one-sided limits as the mean interface position is approached from above and below are denoted respectively by z → 0 + and z → 0 − . The interface drift velocity U 0 and and stress τ 0 at the interface are defined by (1).
The magnitude of the stress at the interface may be expressed in terms of the respective friction velocities v * and u * for the two fluids, which are related by where Let δ 1 and δ 2 , respectively, be intrinsic length scales, such as roughness lengths, that characterize the turbulent boundary layers adjacent to the interface. These friction velocities and characteristic length scales may be used to nondimensionalize the velocities V(z) and U(z) and the height z, giving dimensionless velocitiesV(ζ) andŪ(ζ) that are functions of the dimensionless height ζ, The dimensionless continutiy condition (1) at the interface is then while the dimensionless free-stream velocities arē Consider now the dimensionless velocity deviations in the boundary layers, which satisfy the boundary conditions The problem for the dimensionless velocity deviations depends on the dimensionless free- in the interface condition (10). It will therefore be unchanged if the free-stream velocities V ∞ andŪ ∞ are replaced by a different pair of free-stream velocities that preserve that difference, that is, for modified dimensionless free-stream velocitiesṼ ∞ andŨ ∞ such that The condition (11) will be satisfied if for any choice of ∆V. Now, choose ∆V such that which requires that The dimensionless modified problem with free-stream velocities defined by (12) and (14) has antisymmetric free-stream conditions and, by assumption and scaling, universal turbulent boundary layer structure on both sides of the interface. Given this symmetry and universality, it it consistent to assume that the dimensionless modified total velocity profile is itself antisymmetric across the interface. This in turn implies that the dimensionless modified total velocity vanishes at the interface, that is, The sumŪ (0)+U ∞ in (15) is just the dimensionless interface drift velocityŪ 0 = U 0 /u * for the original problem with dimensionless free-stream velocities V ∞ and U ∞ , and consequently that drift velocity can be computed from (15): The corresponding dimensional interface drift velocity is then: The fundamental universality result for the interface drift velocity is (17). It follows from (17) that the drift velocity relative to the lower fluid is proportional to the velocity of the upper fluid relative to the lower fluid: The ratio of the relative drift speed to the relative free-stream speeds thus depends only on the relative densities of the two fluids ( Fig. 1): It follows from (17)- (19) that and The classical air-water 20 and freshwater-saltwater 19 cases considered below respectively represent these two limits. Intermediate density ratios can be obtained for two-fluid systems consisting of water with other laboratory fluids (Fig. 1).

A. Freshwater-saltwater
It is useful as a reference point to consider briefly the case of two fluids of similar densities, for which the dimensional mean velocity profiles should themselves have a nearly antisymmetric structure about the interface. Experiments on a two-layer, freshwater-saltwater system conducted by Lofquist 19 provide a specific illustration, for which the maximum 7% density difference considered was similar to that for an olive-oil and water system (Fig. 1).
The flow of the lower, salt-water layer was forced by a pump system, and the upper, freshwater layer was driven by momentum transport through the interface. Despite the small density differences between the fluids and the induced turbulence, relatively sharp density interfaces were maintained through the experiments, consistent with the theoretical configuration assumed in Sec. II.
The measured profiles of mean horizontal velocity and its vertical derivative show the anticipated, approximately antisymmetric and symmetric structures, respectively (Fig 2).

B. Air-water
For an air-water laboratory system with limited wave effects, a systematic study of wind-  Wu 20 reports total and "wind-induced" surface currents ranging from less than 10 to over 50 cm s −1 , both of which are in approximate agreement with the general theoretical prediction (20) for surface wind drift (Fig. 3). The "wind-induced" surface current was computed by Wu 20 by removing an estimated surface Stokes' drift velocity from the measured, total surface current. Because of the small fetch and wave-damping baffles in the laboratory apparatus, wave amplitudes were less than 5 cm for the highest free-stream velocities and decreased toward zero for lower free-stream velocities, so the Stokes' drift estimate was relatively small, generally less than 15% of the total surface current. The Wu 20 velocity data show slightly better agreement with the general theoretical prediction (20) than with the alternate empirical relation |U 0 | ≈ 0.53 v * to which Wu 20 compares the laboratory surface drift measurements (Fig. 4), suggesting that the universality prediction (17) may be slightly more robust than the empirical proportionality to the friction velocity.
C. Numerical simulations of Lombardi et al. 26 Lombardi et al. 26  in each case, so the laminar viscous sublayers contained a significant fraction of the shear and were explicitly resolved by the numerical scheme. The kinematic viscosities ν L and ν G for the liquid and gas phases were further taken to be related by ν L = α ν G , where α = (ρ G /ρ L ) 1/2 , for which the resulting R29 case then corresponds approximately to air and water at temperature 320 K and standard atmospheric pressure 26 . The standard boundary conditions (1) were imposed on the instantaneous flow at the interface, which in these simulations was represented as a rigid plane surface and not allowed to deform. Flow was forced in the two layer by equal and opposite dimensional pressure gradients, so that no net momentum was imparted to the system.
For the R1 simulation, with ρ 1 = ρ 2 and α = 1, Lombardi et al. 26 note that the exact antisymmetry (neglecting a small deviation arising from the numerical implementation) of the resulting flow configuration requires that the velocity at the interface be exactly zero, just as in (15). It is clear, then, that for the R1 simulation, the drift velocity relative to the lower fluid is equal to one-half the difference of the free-stream (outer boundary) velocities ( Fig. 1), just as in (19) and (21), and, essentially, as in the freshwater-saltwater example which causes the interface roughly to resemble a no-slip boundary for the gas but a free-slip boundary for the liquid. However, the expression (19) nonetheless provides an accurate first-order prediction of the relative interface drift velocities |U 0 − U ∞ | for both the R10 and R29 simulations (Fig. 1). The error associated with the departure from antisymmetry for these two cases is no more than 2.5/18 ≈ 15%, as the free-stream (outer boundary) velocities (U ∞ , V ∞ ) for these flows are both approximately 18 times the respective friction velocities. For the R29 simulation, this result might be anticipated also from the empirical proportionality |U 0 | ≈ 0.5v * cited by Wu 20 for the air-water case, by which the error might be estimated as 2.5u * /0.5v * ≈ 5α ≈ 15%. Consistent with the Lombardi et al. 26 interpretation, it is shown in Sec. IV B that these small departures from universality for the R29 and R10 simulations can be expressed as corrections to the interface velocity that are logarithmic in the ratio of effective roughness lengths for the outer, turbulent flow regimes.

A. General considerations
The derivation of (17) is general but hides two important assumptions, beyond the standard turbulent universality conditions of homogeneity and stationarity. First, it assumes that the interface geometry is sufficiently simple that averaging can be conducted with respect to a quasi-vertical, curvilinear-coordinate distance from the interface. Second, it relies on a dimensional analysis argument that is incomplete if either or both boundary layer flows are characterized, respectively, by more than a single intrinsic length scale. Some level of violation of these assumptions is inevitable in every observable fluid system, and the exact results (17)- (19) are thus best regarded as approximations or a priori estimates.
The assumption on the interface geometry should be satisfied as long as the wave activity at the interface is sufficiently weak that the fluctuating interface height has the form of a continuous, single-valued function of horizontal position and time. It is reasonable to expect that it could more generally be satisfied as long as the interface remains a continuous surface, that is, prior to energetic wave-breaking events that cause either bubble and spray formation, or irreversible mixing between the layers, depending on the fluid physics. Whether the result may provide useful guidance after the onset of energetic wave-breaking could be assessed empirically and presumably would depend on the fraction of the surface affected by the wave-breaking events. When the interfacial shear is sufficiently large relative to the density difference, frequent wave breaking may destroy the interface, in which case the definition of interface velocity itself becomes difficult or impossible. It is important to note, however, that a role of interfacial waves in supporting the stress τ 0 at the interface is not excluded, as the result (17) does not depend on any direct assumptions regarding the mechanisms of momentum transport.
An additional intrinsic length scale that would violate the dimensional analysis assumption will in general arise whenever the outer boundary layer scales are constrained independently of δ u and δ v . In laboratory flows, for example, the outer scales may be controlled by the dimensions of the apparatus. In rotating flows, with rotation rate measured by Ω, the outer scales may be constrained by the inertial depths (δ v , δ u ) = (v * , u * )/Ω, for which which fit the R29 numerical profiles when ζ G 0 = 0.14 and ζ L 0 = 0.42 (Fig. 5). It follows that and therefore, withŪ ∞ = −V ∞ for the Lombardi et al. 26 simulations, For the fitted roughness lengths (ζ G 0 , ζ L 0 ) = (0.14, 0.42) and α = 1/29.9464 for the R29 simulation, the correction (24) givesŪ 0 = −2.65. This is close to the R29 numerical result Table II of Lombardi et al. 26 , as it must be, given the fitting procedure.
for the linear region radze logarithmic law e also shown. The proves and, with good ap-CH~see Table I!. The ences in boundary con-, continuity of velocity '' rofile on the liquid side ose to the interface the that in the wall region tly lower than the CH a significantly lower acial velocity. The dife to differences in the olds stress (u8w8) proat for the case CH, as alue of the interfacial difference in the mean profiles were identical to be for R51! then, set up, the interfacial alues of the interfacial nt on the two sides of nsionalization is done , from the definition of e II shows how the inuld ideally be zero for inaccuracy of about s of the mean streamthe reference case, CH, profile is reported with e. not subtracted out as profiles are very simi-H. At z G

1
.170, about the liquid side, on the nd R10 runs are differterfacial velocity is left r in the log-law region 0. However, there is an h useful in understand-~b! shows that the vevery similar to the case acted out, whereas From the view point of the liquid, the surface parallel fluctuations can proceed relatively unimpeded as the gas has low inertia. On the other hand, the gas side fluctuations are impeded by the high inertia of the liquid. Thus the gas sees the interface as something like a no-slip boundary, whereas for the liquid it is similar to a free-slip boundary. Consequently, interface-parallel fluctuations on the liquid side can be larger than in wall turbulence. Due to the significantly larger streamwise velocity fluctuations the Reynolds stress increases and the velocity profile is therefore altered in the buffer region. The effect of the boundary dies away in the log-law region.

1649
Lombardi, De Angelis, and Banerjee rameter R, defined as the square root of the density ratio between the fluids~if applicable!.
The run CH is a reference case where a ''no-slip'' boundary condition is applied at the interface, i.e., the fluids are uncoupled and the flow corresponds to that in a channel with one slip wall and one no-slip wall.
R is an indicator of the degree of dynamical coupling between the phases, as evident from~11!. Previous studies 12-16 assumed a high value of the density ratio between the fluids and uncoupled the flows, i.e. they analyzed the limiting case,
From these runs, it is possible to determine the main features of turbulence in the vicinity of a sheared fluid-fluid interface, and how the structures vary as the density ratio varies. Intuitively, it may be expected that the interactions increase as the density ratio approaches unity.
The Reynolds number, Re, has been kept at 60.4 for both subdomains~see Table I!. This is the same value used by Lam and Banerjee 16 These are the same values chosen by Lam and Banerjee, 16 for this Reynolds number, in the streamwise and spanwise directions.
To satisfy the constraints on h and h 1 discussed above, the two fluids must be such as to fullfill the condition This is close to many physical situations. For example, the R29 case~with R529.9464! corresponds to air-water at atmospheric pressure and 320 K roughly, of interest for transport processes in environmental systems. The higher density ratios may occur in industrial equipment, e.g. stratified flow in oil-gas pipelines. There is, of course, also intrinsic interest in investigating the effect of R on interfacial turbulence characteristics in limiting cases. R1 represents the minimum in the range of R, with equal density fluids, separated by a plane interface across which shear stresses can be transmitted, but which is not allowed to deform. This is not a physically realizable situation but of interest for purposes of comparison. Run R10 corresponds to an intermediate condition.
Considering~18! and~21!, the interfacial shear stress boundary conditions~12! simplify to which, of course, holds only in the particular case considered here of flow with the same surface-normal extent of the subdomains in shear-based units.
In setting up the problem, we have adopted the convention that the top fluid is lighter than the bottom, though the convention may, just as well, be reversed, since the interface is not allowed to deform in the present study.

A. Mean velocity profiles
Figures 5~a! and 5~b! show the profile of the computed mean velocity~streamwise component! for case R29~normalized by the interfacial shear velocity, u t ), the plotted quantity being: ū 1 5ū/u t . As mentioned earlier, this case is similar to air-water flow at atmospheric conditions. The interface is at z 1 50. The non-dimensional value of the mean interfacial velocity has been subtracted out from the profile, to allow a comparison with the reference run, CH, where the interface corresponds to a ''no-slip'' wall.~As discussed later, the interfacial velocity is not exactly zero because, in spite of the same non-dimensional pressure gradient -in opposite directions-in each phase, the flowrates are somewhat different, due to the Reynolds stresses on each side being Velocity profiles for the R10 numerical simulations (Fig. 5) are nearly identical to those for the R29 simulations, but according to 95. This relatively small inferred change in effective roughness lengths is consistent with the near-equality of the R10 and R29 velocity profiles in the logarithmic outer layer (Fig. 5).

C. Corrections for the Wu 20 air-water flow
It might be assumed from the agreement (Fig. 3) of the surface drift measured by Wu 20 with the theoretical prediction (20) that the symmetry and universality assumptions are satisfied nearly exactly in this laboratory flow. The validity of this inference can be explored by examining solutions of an explicit analytical model of the laboratory flow, which reveals that the apparent close agreement masks two compensating, small but non-negligible, departures from universality, each of order 15% in α. These evidently arise from wind-wave disequilibrium and from differing air and water ratios of outer boundary layer scales to effective interfacial roughness lengths. The former represents a departure from homogeneity that effectively modifies the stress condition at the interface, while the latter results in a correction that is logarithmic in the additional dimensionless parameter. In each fluid, the pressure gradient force is balanced by the divergence of a turbulent stress that is parameterized by an eddy viscosity, where the latter is taken to increase linearly with distance from the interface and the top and bottom of the tank sufficiently near to those boundaries, reaching a constant maximum sufficiently far from those boundaries. For air or water roughness lengths sufficiently small that the corresponding eddy viscosity κv * z a 0 or κu * z w 0 is less than the air or water molecular viscosity ν a or ν w , a viscous sublayer region is included adjacent to the interface, as in the analytical model of the Lombardi et al. 26 numerical simulations (App. B,C).
The dimensional model solutions have an air velocity maximum near 0.20 m height, consistent with the air velocity profiles shown by Wu 28 , and a water velocity with narrow top and bottom boundary layers and a weakly sheared interior passing through zero near mid-depth (Fig. 6). The boundary layer widths in both air and water are of order 0.1 m.
The velocity difference across the upper water boundary layer, adjacent to the interface, is several times greater than that across the lower water boundary layer.
The dimensionless model surface-relative velocity profiles, [V (z/z a 0 ) − U 0 ]/v * for the air and [U (z/z w 0 ) − U 0 ]/u * for the water, are effectively antisymmetric up to dimensional distances of order 10 3 roughness lengths from the interface, and approximately antisymmetric for distances approaching the height of the free-stream air velocity maximum (Fig. 7). This universal and antisymmetric structure is consistent with the assumptions from which the general theoretical results (17)-(20) are derived. As anticipated, the universal structure exhibited by these solutions is a classical logarithmic wall boundary layer (Fig. 7), consistent with the measured laboratory velocity profiles on both sides of the interface 20,28 .
The logarithmic boundary layer structure and the theoretical result (20) can be combined to yield a prediction for the empirical proportionality |U 0 | ∝ v * reported by Wu 20 and previous authors. The universal boundary layer solution takes the form |V|/v * ≈ (1/κ) ln(z/z a 0 ). If the free-stream velocity V ∞ is approximated by the boundary layer solution at a height z a ∞ , with U ∞ ≈ 0 from the mean or mid-point water velocity, then (20) becomes with direction determined by the free-stream velocity. Madsen 30 , (his eq. 40, Sec. 4), used similar but converse reasoning to obtain a 3% proportionality of wind drift to 10-m wind speed from his model prediction of wind-drift dependence on friction velocity for the oceanatmosphere case. The constant of proportionality in (25), (α/κ) ln(z a ∞ /z a 0 ), depends only on the ratio of the free-stream velocity height to the roughness length, and for the oceanatmosphere case is directly related to the standard 10-m neutral drag coefficient 31 .
For the model solutions, z a ∞ /z a 0 ≈ 10 3 (Fig. 7), so that (25) gives |U 0 | ≈ 0.60 v * , similar to the relation |U 0 | ≈ 0.53 v * reported by Wu 20 . Direct evaluation of the ratio |U 0 |/v * from the model solutions also gives values that roughly match the individual observed values at each free-stream velocity, which range from 0.4 to 0.7 (Fig. 4). The larger variation with V ∞ (Fig. 4a) of the proportionality constant in the relation U 0 ∝ v * , compared to that for the equivalent α in (20), evidently arises from the dependence seen in (25) on the roughness length z a 0 , which in turn depends on V ∞ (Fig. 8). values, then the model can be seen to have a unique dimensionless solution that scales with the pressure-gradient forcing, and the ratios U 0 /V ∞ and U 0 /v * are both constant, as can be verified by numerical solution.
One departure from universality is indicated by the difference in the Wu 20 estimates of interfacial stress obtained from the air and from the water velocity profiles, for which the water friction velocities are systematically smaller than would be predicted by (3) from the observed air friction velocities (Fig. 9). As noted by Wu 20 , this difference likely indicates that a portion of the stress from the air is absorbed directly by the wave field, which was not in local equilibrium at the mid-tank measurement point, but developed downstream in response to the surface forcing until damped by the downstream baffles.
A second departure is indicated by the difference in the modeled and measured water roughness lengths: for this set of solutions, the former are one to two orders of magnitude larger than the latter (Fig. 8b). With other parameters held fixed, variations in the model   water roughness length z w 0 have a much larger relative effect on the model surface drift velocity U 0 than on the air free-stream velocity V ∞ and friction velocity v * (Fig. 10). For z w 0 > ν w /(κu * ), the drift velocity decreases linearly with the logarithm of z w 0 ; for smaller z w 0 , the dependence is weaker, because the viscous sublayer effectively prevents the surface from becoming arbitrarily slippery as z w 0 → 0. Note that, perhaps surprisingly, the measured water roughness lengths are smaller than the measured air roughness lengths, rather than larger as in the numerical simulations of Lombardi et al. 26 .
For the cases with V ∞ ≤ 9 m s −1 , the model fit can be improved by changing the stress condition (2) to reflect an assumption that roughly one-quarter of the stress goes into the wave field rather than the mean current, consistent with the estimates by Wu 20 , and simultaneously reducing the water roughness lengths to z w 0 = 0.1 z w 0 (Figs. 9,8b). These adjustments bring the model and observed water friction velocities and roughness lengths into better agreement, while maintaining the model fit to the observed V ∞ , v * , and U 0 .
The modified stress condition is imposed in the model by setting the value of α to α = 0.85 α ≈ 0.029: because the fluid densities appear only in α, which enters only through (2), this is equivalent to increasing the water density by 40%, so that the stress results in a reduced water velocity. For V ∞ > 9 m s −1 , the water roughness lengths are reduced only to z w 0 = 0.33 z w 0 in these modified solutions, as a greater reduction results in systematic over-prediction of the drift velocity, for reasons that have not been determined but may be related to wave effects, either on the dynamics or on measurement error.
An analytical expression for the corresponding leading-order correction to the predicted with dimensionless equivalents where ζ = z/z a 0 for z ≥ 0 and ζ = z/z w 0 for z ≤ 0. Let z a ∞ and z w ∞ be levels at which V and U approach their respective free-stream values, V ∞ and U ∞ , with ζ a ∞ 1 and ζ w ∞ 1 the corresponding dimensionless levels. Then and from which it follows that If |ζ w ∞ | = ζ a ∞ , so that the symmetry assumption holds, (30) reproduces (18). If, however, the water roughness length is replaced by a smaller value z w 0 = βz w 0 , then |ζ w ∞ | is replaced by |ζ w ∞ | = |ζ w ∞ |/β = ζ a ∞ /β and (30) yields the leading-order correction to (18), With u * held constant at the value obtained for the solution with |ζ w ∞ | = ζ a ∞ , consistent with weak dependence of friction velocity on z w 0 (Fig. 10b), the approximation (31) quantitatively reproduces the leading-order linear dependence of the drift velocity on the logarithm of the water roughness length that was found in the numerical solutions (Fig. 10). The dependence is such that the drift velocity increases with decreasing water roughness length, and decreases with increasing water roughness length. Note that the stress is constant in the high-shear near-surface log layers, and consequently the surface stress remains constant while the drift velocity changes in response to the variations in roughness length.
a modified form of (20). The improved prediction (32) for the wind drift, which includes the leading-order correction for asymmetry of the dimensionless velocity profiles, may be compared with (20). For the model solutions with α replaced by α = 0.85α and z w 0 replaced by z w 0 = 0.1z w 0 (Figs. 8,11), β = 0.1 and ζ a ∞ ≈ 10 3 , so that (32) givesα ≈ 1.3 α ≈ 1.1α. Thus, for these solutions, the effects of the dimensionless asymmetry and of the partial absorption of the stress by the wave field essentially compensate, giving a predicted effective valueα for the proportionality of wind drift to free-stream velocity that happens nearly to coincide with the true value of α. Bothα and α then appear to predict accurately the measured wind drift, especially for V ∞ < 9 m s −1 (Fig. 3).

V. SUMMARY
The theoretical argument leading to the results (17)- (21) for the drift velocity on the interface relies on a classical assumption of turbulent universality and Reynolds' number similarity, which asserts that the statistically steady and horizontally homogeneous boundary layers on either side of the interface must be controlled by the same universal turbulent processes. This equivalence then implies that the mean velocity deviation profiles, when scaled by the respective friction velocities and intrinsic length scales, will be antisymmetric across the interface. This condition is sufficient to determine the interface drift velocity in terms of the free-stream velocities. The resulting expression (17) provides a simple, physically intuitive, and quantatively useful first-order prediction of the interface drift velocity for two-fluid systems with general density ratios. The comparisons with classical laboratory experiments on freshwater-saltwater and air-water systems and with previous numerical simulations of coupled gas-liquid flow illustrate this result. In addition, they provide insight into the amplitude and physical origins of corrections to the predicted drift velocity that will typically result from departures from universality, which are to be expected in all realizable fluid systems.
Among the many possible combinations of fluids to which the result (17) may be relevant, the air-water system is of particular interest, especially in the context of surface wind drift in the ocean-atmosphere system. The classical empirical rule for ocean wind drift reported already by Keulegan 16 and Van Dorn 17 has remained remarkably consistent over many decades and in a wide range of environmental conditions: Wu 20 states that it is "commonly accepted" that the total wind drift is about 3% of the wind velocity at long fetches; Madsen 30 and Weber 32 cite a "commonly employed rule of thumb" that oil slicks are advected with a velocity which is 3% of the wind speed, in a direction that is approximately 10 • or 15 • to the right of the wind; Morey et al. 33 characterize the wind drift as a correction to numerical model surface velocity that is "typically 3%" of the wind speed and directed "to the right of the wind direction;" while for the contemporary operational ocean model described by Zelenke et al. 18 , the suggested proportionality is "typically about 3% of the wind speed," but may range from 1% to 4%, depending on conditions. This long-standing empirical rule nonetheless remains without definitive observational confirmation or accepted theoretical explanation.
It is striking that the universality predictions (17)- (20) essentially reproduce this empirical rule, and consequently it is also tempting to hypothesize that the approximate universality of the turbulent dynamics gives a plausible explanation for its stability and reliability.
In this interpretation, the apparent rotation to the right of the surface (10-m) wind in the northern hemisphere would be understood as arising from the rotation of the surface (10-m) winds to the left of the geostrophic wind, the rotating boundary layer analog of the free-stream velocity 34 , which by (21) would determine the direction of drift velocity.
Arguing against this interpretation, especially, is the ubiquitous role of wave-breaking on ocean surface dynamics at all but the lowest wind speeds, and the expected strong effect of surface waves on near-surface Lagrangian motion, with much recent work on the ocean surface drift problem focused instead on estimating a wave-driven Stokes' drift at or near the surface [11][12][13][14][15]33 .
In addition to its long-standing importance for such practical problems as pollutant dispersal, search and rescue, and navigation, the determination of surface ocean drift has recently drawn renewed scientific attention. Simultaneous measurement of ocean surface winds and currents has been identified, for example, as a priority for future satellite observations by the 2018 Decadal Review 35 , which cites their importance in a variety of Earth science contexts, including the determination of air-sea momentum exchange, ocean upwelling, upper ocean mixing, and sea-ice drift. An airborne Doppler scatterometer has been developed that simultaneously measures surface stress and currents in the upper few centimeters of the ocean 36 , and a design for a satellite version of this instrument exists 37 . Novel in-situ technologies are also in development that promise better direct measurements of currents within centimeters of the ocean surface 2,33 . Improved understanding of the dynamics at the atmosphere-ocean interface is sure to follow from these developments, and from improved theoretical insight into the controlling processes. It is hoped that the present contribution may prove to be useful in this context and help to stimulate further work on these topics.      For interface air roughness length z a 0 ≥ ν a /(κv * ), where ν a = 1.5 × 10 −5 m 2 s −1 is the kinematic viscosity of air at 20 • C, the model equations for the air (z > 0) are: where τ a is the turbulent stress, G is the pressure gradient, ρ a is the air density, and A V (z) is the air eddy viscosity, defined by: In (B2), κ = 0.4 is the von Kármán constant, v * and v * H are the wind friction velocities at the interface and the tank top, z a 0 and z H 0 are roughness lengths associated with the interface and the tank top, z a 1 and z a 2 are constants, and H is the height of the tank top. For interface water roughness length z w 0 ≥ ν w /(κu * ), where ν w = 1.0 × 10 −6 m 2 s −1 is the kinematic viscosity of water at 20 • C, the model equations for the water (z < 0) are analogous: where τ w is the turbulent stress, F is the pressure gradient, ρ w is the water density, and A U (z) is the water eddy viscosity, defined by: In (B4), u * and u * D are the water friction velocities at the interface and the tank bottom, z w 0 and z D 0 are roughness lengths associated with the interface and the tank bottom, z w 1 and z w 2 are constants, and D is the depth of the tank bottom. In order that the eddy viscosities be continuous, the constants z a 2 and z w 2 are computed from the conditions The boundary conditions are that the velocities vanish at the top and bottom of the tank, supplemented by continuity of the velocities and the stress at the interior fluid points z a 1 , z a 2 , z w 1 , z w 2 and at the interface z = 0. In addition, the condition of no depth-integrated water flow is imposed: The momentum equations (B1) and (B3), with (B2) and (B4), can be integrated analytically in each subregion, subject to unknown constants to be determined by the boundary and continuity conditions. For the solutions considered here, imposed air pressure gradients G < 0 are chosen, so that V ≥ 0, dV /dz < 0 at z = H, dV /dz > 0 and dU/dz > 0 at z = 0, and dU/dz < 0 at z = −D. The friction velocities, which are determined as part of the solution, then satisfy Vertical integrals of the momentum equations (B1) and (B3) yield while (B8) implies (4) with ρ 1 = ρ a , ρ 2 = ρ w , so that v * H , u * , and F can be computed directly from G, v * , and u * D . Solutions are obtained by specifying the air pressure gradient G and the roughness lengths at the interface and the top and bottom boundaries, and solving simultaneously, by numerical iteration, for the the water pressure gradient F and the friction velocities.
The model equations may be integrated to find explicit solutions in each region, which satisfy the continuity conditions on velocity and stress at z = {z a 1 , z w 1 }, the continuity conditions on stress at z = 0, and the no-slip boundary conditions at z = {−D, H}, as follows: z a 1 ≤ z ≤ z a 2 : V (z) = V (z a 1 ) + 1 κv * (z a 0 + z a 1 ) For interface roughness lengths z a 0 < ν a /(κv * ) or z w 0 < ν w /(κu * ), a corresponding laminar viscous sublayer analogous to that in (C1) is included adjacent to the interface, with continuity of velocity and stress imposed at the transition points z = {z a lam , −z w lam }, where z a,w lam = {ln[ν a /(κv * z a 0 )] ν a /(κv * ), ln[ν w /(κu * z w 0 )] ν w /(κu * )}. The log-layer solutions in the regions z a lam < z < z a 1 and z w 1 < z < −z w lam are then also modified from (B11) and (B14) as in (C1). Note that ν w /ν a = 0.067 ≈ 2α, so these solutions do not satisify the additional symmetry condition ν L = αν G imposed by Lombardi et al. 26 .
By (B10), v * H , u * , and F can be computed directly from G, v * , and u * D , using the continuity of stress at z = 0 and vertical integrals of the momentum equations. The velocity profiles and roughness lengths at the tank top and bottom were not measured by Wu 20,28 , so the corresponding model roughness lengths were likewise chosen simply to be proportional to the respective interface roughness lengths, with z H 0 = 0.1 z a 0 and z D 0 = 0.2 z D 0 . For given G, an iterative numerical method may be used to find the values of v * and u * D for which the solutions (B11)-(B16) are continuous at z a 2 and z w 2 and the integral condition (B6) is satisfied. The method was implemented by computing three independent estimates of the unknown drift velocity U 0 in (B11)-(B16), from the integral relation (B6) and the continuity conditions at z a 2 and at z w 2 , and then minimizing the squares of the differences of two pairs of these estimates.
An initial set of solutions approximately matching the Wu 20 observations of free-stream velocity and air friction velocity were obtained by iteration, primarily through adjustment of the air pressure gradients G and air roughness lengths z a 0 , with the first guesses for the air roughness lengths taken from the values reported by Wu 20 . For these solutions, the water roughness lengths were set according to z w 0 = 4 z a 0 for V ∞ ≤ 9 m s −1 and z w 0 = 3 z a 0 for V ∞ > 9 m s −1 . Given the experimental uncertainties, the iteration was conducted manually and halted when approximate agreement of the model free-stream and air friction velocities with the measurements was achieved.
Appendix C: Analytical model solutions for Lombardi et al. 26 flow For the analytical model of the Lombardi et al. 26 flow, the domain is divided into a laminar sublayer region 0 < |ζ| ≤ ζ lam with unit dimensionless molecular viscosity and a turbulent outer region with dimensionless eddy viscosity κ(ζ κ + ζ). The transition point ζ lam and the roughness-length parameter ζ κ are related by the requirement that the eddy and molecular viscosities be equal at the transition point, where continuity conditions are imposed on the velocity and stress; i.e., that κ(ζ κ + ζ lam ) = 1. The resulting dimensionless model velocity profilesV andŪ in the gas (G) and liquid (L) arē where ζ = zv * /ν G = zu * /ν L (for ν L = αν G , as assumed), ζ 1 = 2 √ 2 Re = 170.8, andŪ 0 is the dimensionless interface velocity. The outer boundary conditions (dV /dζ = dU/dζ = 0 at ζ = ζ 1 ) cited by Lombardi et al. 26 would suggest F = 1/ζ 1 but a value of F = 0.1/ζ 1 was used to obtain a better fit to the numerical velocity profiles in Figure 5. For the log-layer only model, these equations are replaced bȳ for which the dimensionless velocities relative to the interface,V =V −αŪ 0 andŪ =Ū −Ū 0 , are purely logarithmic with effective dimensionless roughness lengths ζ G 0 and ζ L 0 . The liquidside profiles are reflected about ζ = 0 in Figure 5, i.e., |Ū (ζ)| and |Ū (ζ)| are shown vs. |ζ|.