Articles | Volume 19, issue 1
Research article
24 Jan 2023
Research article |  | 24 Jan 2023

Extension of Ekman (1905) wind-driven transport theory to the β plane

Nathan Paldor and Lazar Friedland

The seminal Ekman (1905) f-plane theory of wind-driven transport at the ocean surface is extended to the β plane by substituting the pseudo-angular momentum for the zonal velocity in the Lagrangian equation. When the β term is added, the equations become nonlinear, which greatly complicates the analysis. Though rotation relates the momentum equations in the zonal and the meridional directions, the transformation to pseudo-angular momentum greatly simplifies the longitudinal dynamics, which yields a clear description of the meridional dynamics in terms of a slow drift compounded by fast oscillations; this can then be applied to describe the motion in the zonal direction. Both analytical expressions and numerical calculations highlight the critical role of the Equator in determining the trajectories of water columns forced by eastward-directed (in the Northern Hemisphere) wind stress even when the water columns are initiated far from the Equator. Our results demonstrate that the averaged motion in the zonal direction depends on the amplitude of the meridional oscillations and is independent of the direction of the wind stress. The zonal drift is determined by a balance between the initial conditions and the magnitude of the wind stress, so it can be as large as the mean meridional motion; i.e., the averaged flow direction is not necessarily perpendicular to the wind direction.

1 Introduction

The seminal theory of wind-driven transport at the ocean surface was developed about 120 years ago by the Swedish oceanographer Vagn Walfrid Ekman for the highly idealized case of constant Coriolis frequency – the f plane. The Ekman (1905) theory addresses the downward-spiraling horizontal velocity in the ocean's surface and its vertical integral – the transport. Ekman's elegant solution of the problem has become textbook material in physical oceanography, dynamical meteorology, and geophysical fluid dynamics (see, e.g., Gill1982; Pedlosky1987; Vallis2017). For uniform wind stress the dynamics on the f plane consist of two parts: a steady flow to the right and left of the wind direction in the Northern and Southern Hemisphere and inertial oscillations (of frequency f0 – the constant Coriolis frequency). However, though it is one of the cornerstones of atmosphere and ocean dynamics, the theory was never extended to include the latitudinal increase in the Coriolis frequency, known as the β effect, which is the focus of the present study. In contrast to the β plane, in spherical coordinates the theory of wind-driven transport was studied numerically in Constantin and Johnson (2019) and Paldor (2002), but due to the complexity of the governing equations in these coordinates, the numerical solutions have not yielded analytic understanding. With the wind-driven dynamics on the f plane fully understood and quantified, the β plane offers an in-between set-up wherein analytical insight can complement the numerical solutions.

For given wind stress forcing, the known general differences between the dynamics on the f plane and β plane heuristically suggest that the extension of Ekman's transport theory to the β plane should include the following qualitative elements.

  1. An increase or decrease in mean meridional velocity for an eastward- or westward-directed stress due to the decrease or increase in Coriolis frequency when the water column moves southward or northward.

  2. The frequency of oscillation about the mean velocity should decrease or increase (so oscillation period should increase or decrease) due to the decrease or increase in Coriolis frequency along the trajectory (for an eastward-directed stress in the Northern Hemisphere, while the opposite changes occur for westward-directed stress and in the Southern Hemisphere).

  3. Since the oscillation's frequency and amplitude are inversely correlated (energy flux is unchanged) a decrease in frequency should lead to an increase in amplitude and vise versa.

  4. Since inertial oscillations that form a perfectly circular motion on the f plane drift westward on the β plane, the averaged zonal motion should drift to the west. A heuristic reasoning of the westward drift in terms of the change in the radius of the inertia circle was proposed by Von Arx (1964), and complete quantitative theories of the drift were developed in Ripa (1997) and Paldor (2007).

The numerical solutions of the governing Lagrangian equations (see Sect. 2 below) shown in Fig. 1 fully confirm the first three expectations listed above but contradict the fourth one – for both westward (right panel) and eastward (left panel) stresses, the trajectories drift to the east. From the particular example shown in Fig. 1 it is unclear whether the eastward transition is a general feature of the wind-driven dynamics on the β plane or a specific occurrence related to the particular choice of initial conditions and/or parameter values.

Figure 1The (longitude, latitude) trajectories of water columns at the ocean surface subject to westward-directed (b) and eastward-directed (a) wind stress on the f plane (blue curves) and on the β plane (red curves). The time unit is the inverse of the mean Coriolis frequency, and the longitude and latitude distances are scaled on Earth's radius. The value of b (scaled β) corresponds to 30 latitude. The scaling of the wind stress (τx) is detailed in Sect. 2. Both trajectories start from (x,y)=(0,0) located at the bottom right point in the right panel and at the upper left point in the left panel.


In addition to resolving the issue of the zonal drift and quantifying the various rates of changes, the present study also addresses the equatorial problem that exists only on the β plane. This equatorial issue can be described as follows: an eastward-directed stress in the Northern Hemisphere forces a net southward-directed mean flow, which, on the β plane, is accompanied by a decrease in the Coriolis frequency. Thus, at some time the wind-forced water column must find itself at a latitude at which the Coriolis frequency vanishes – the Equator. From that point onward the water column is subject to non-rotating dynamics and must move eastward at an accelerated velocity. In the rest of this work we will estimate the time it takes the water column to change its dynamics from rotating to non-rotating and analyze how the two dynamical regimes connect with one another.

The work is organized as follows: in Sect. 2 we nondimensionalize the governing Lagrangian equations and simplify them by substituting the pseudo-angular momentum for the zonal velocity. The simplified system is analyzed in Sect. 3, and the work concludes with a discussion and summary in Sect. 4.

2 The nondimensional model

The time-dependent trajectory of a column of water in the surface Ekman layer forced by the overlying uniform wind stress on the f plane is a fundamental problem of physical oceanography that is fully described in most textbooks (Gill1982; Pedlosky1987; Vallis2017). The governing Lagrangian equations that describe the dynamics of vertically averaged horizontal velocity components consist of the momentum equations in the zonal and meridional directions and the (trivial) relations between these velocity components as well as the changes in the coordinate of the moving column, i.e.,

(1) d x d t = U , d y d t = V , d U d t = f V + τ x H ρ , d V d t = - f U .

Here τx is the uniform zonally directed wind stress (which is positive or negative for eastward- or westward-directed wind, respectively), ρ is the water density, H is the depth (thickness) of the layer, f=f0+βy is the Coriolis parameter (where f0=2Ωsin(ϕ0), β=2Ωcos(ϕ0) /Re with Re and Ω – Earth's radius and rotation frequency, respectively, and ϕ0 – the latitude at which the plane is tangential to Earth), U and V are the vertically averaged horizontal velocity components in the eastward and northward directions, respectively, and x and y are the respective coordinates in these directions. The only added complication of this system relative to that studied in detail in, e.g., Ch. 9 of Gill (1982), is that here the Coriolis frequency, f, in the momentum equations is y-dependent.

The four-dimensional system (Eq. 1) can be easily integrated numerically, but the general properties of its solutions can be best deciphered by reducing the number of its free parameters. This is done by scaling time, t, on 1f0 and x and y on Re so the velocity scale is f0Re. With this scaling the nondimensional Coriolis frequency is 1+by, where b=βRef0=cot(ϕ0) is the nondimensional β. The system is further simplified by replacing U by the pseudo-angular momentum, defined as D=U-y(1+b2y) in nondimensional units. As was shown by Paldor (2007), when τx=0, i.e., in the inertial case, D is conserved. We note that in spherical coordinates the conservation of angular momentum, which is the spherical counterpart of D, yields a simple relation between the zonal velocity and the latitude (Paldor2001). Formally, a similar quantity relating the zonal velocity, U, and the meridional coordinate, y, can also be derived in Cartesian coordinates, but, unlike spherical coordinates, this conserved quantity is not the angular momentum. With these changes the system (Eq. 1) transforms to


Here t, x, y, and V denote the nondimensional counterparts of the dimensional variables denoted by the same symbols in Eq. (1), and, as explained above, D=U-y(1+b2y) is the nondimensional pseudo-angular momentum. Equation (4) confirms that D is indeed conserved when Γ=0. The solutions of this system are determined by the four required initial conditions and the two parameters b=βRef0= cot(ϕ0), the nondimensional β, and Γ=τxHρf02Re, the constant nondimensional surface wind stress. The value of b at ϕ0=30 is 1.75. For realistic values (τx/ρ2×10-4 m2 s−2, f0=10-4 s−1, and H=30 m), the value of Γ=10-3, so the theory should be applicable to b of O(1) and Γ≪1. The sign of Γ is that of τx – positive for eastward-directed stress and negative for westward-directed stress.

We solve this system by starting at the origin of the β plane, i.e., x(0)=y(0)=0, and assume that the initial V(0) and D(0)=U(0) are sufficiently small. The numerical solutions presented below are initiated with D(0)=0 and V(0)≠0. However, the definition of D implies that trajectories emanating from D(0)≠0 and y(0)=0 can also be calculated starting from D(0)=0 and a suitable y(0)≠0. Note that the choice D(0)=0 does not restrict the generality of our solutions since the shift of time from t to t=t+D(0)/Γ yields D(t=0)=0, so D(0)=0 can be assumed. The analysis of the solutions of Eqs. (2)–(5), including numerical examples, are presented in the next section.

Figure 2A schematic demonstration of the change in the potential Φ(y,t) for b=0.1 and Γ=0.1 at t=0,10,20,,100. The direction of increase in time is indicated by the black arrows for t<tcr (a) and t>tcr (b). The minima, ym+, of the potentials are indicted by red circles.


3 Analysis

The analysis of Eqs. (2)–(5) begins with the (V, y) subsystem, i.e., Eqs. (3) and (5) along with the (trivial) solution Dt of Eq. (4). The derived solution of y(t) will then be substituted in Eq. (2) to yield the zonal propagation speed. First, we combine Eqs. (3) and (5) to the single second-order equation

(6) d 2 y d t 2 = - ( 1 + b y ) D + y 1 + b 2 y .

We will discuss solutions of this equation for initial conditions in the vicinity of y=dy/dt=V=0 and assume that Γ is sufficiently small (for the smallness condition see Eq. A3 in the Appendix). We proceed by rewriting Eq. (6) as

(7) d 2 y d t 2 = - Φ ( y , t ) y ,


(8) Φ ( y , t ) = 1 2 Γ t + y 1 + 1 2 b y 2 .

Equation (7) describes the dynamics of a quasi-particle in a slowly (for small Γ) time-varying quasi-potential well Φ(y,t). In Fig. 2 we illustrate this potential for Γ=b=0.1 at times t=0,10,20,,100. The minima of these potentials, denoted collectively by ym, are given by the three roots of

(9) Φ y | y = y m = Γ t + y m 1 + 1 2 b y m ( 1 + b y m ) = 0 .

Two cases should be considered depending on time being below or above the critical time

(10) t cr = 1 2 b Γ .

For t<tcr, there are two minima defined by Γt+ym1+12bym=0, i.e.,

(11) y m ± = 1 b ( - 1 ± 1 - 2 b Γ t ) ,

while for t>tcr, there is a single minimum located at

(12) y m 0 = - 1 b .

Figure 3 shows the numerical solutions of y(t) when the column originates near ym+; i.e., x=y=D=0 and V=0.002. As predicted, the exact solution of y(t) oscillates with small amplitude (that increases with the value of V(t=0)) about the evolution curves of ym+ and ym0 shown by the black curves. The direction of evolution of ym+(t) and its transformation to the constant ym0(t) at t=tcr correspond to the evolution of the red circles in Fig. 2 for t<tcr (a) and t>tcr (b). The averaged numerical solution is expected to deviate appreciably from the simple scenario shown here only for high oscillation amplitude near ym+.

The main idea of the following analysis is that since the system starts near y=dy/dt=V=0, i.e., near the minimum of the potential, ym+, by the adiabatic theory (see pages 531–535 in Goldstein1980), it will stay near this minimum for t<tcr. At t=tcr, ym+ transforms into ym0, and therefore at all t>tcr, the system remains near ym0. Thus, the column remains near the minimum of Φ at all times, while this minimum slowly decreases for t<tcr and stays constant for t>tcr. Since the trajectory originates near the minimum ym+ and since for small Γ the variation of the potential is slow (see Eq. 9), we expect the solution for y to be of the form

(13) y = y m ( t ) + δ y ,

where ym(t) starts at ym+ and later (i.e., at t=tcr) transforms into ym0, and δy is a small perturbation. We substitute this form of solution into Eq. (6) and rewrite the resulting equation as

(14) d 2 δ y d t 2 = F - ω 0 2 ( t ) δ y - A δ y 2 - B δ y 3 ,

where F=-d2ym/dt2 is an inhomogeneous forcing term, and the coefficients of the other three terms on the right-hand side (RHS) of this equation are




In the present model, Eq. (11) implies that for t<tcr, d2ym+/dt2=-bΓ2(1-2bΓt)-3/2, so according to Eq. (15) F=-d2ym+/dt2=bΓ2/ω03>0. The second term on the RHS of Eq. (14) describes linear oscillations having slowly varying frequency ω0(t), while the third and fourth terms represent the effect of small anharmonicity of the potential well near the minimum. Note that for ym=ym+ the term ym+(t) in Eq. (13) describes slow monotonic variation of the latitude shown by the black arrows in Fig. 2 at t<tcr. No such variation exists at t>tcr since then ym=ym0= const. As will be shown below, the nonlinear terms in Eq. (14) mostly affect the zonal drift in x.

Figure 3Numerical solutions of y(t) in Eqs. (2)–(5) starting from x=0=y=D and V=0.002. The black curve shows the evolution of ym+(t) given by Eq. (11) for t<tcr and of ym0(t) given by Eq. (12) for t>tcr. The values of b=2 and Γ=0.005 used here imply that the change between the two approximate solutions occurs at tcr=12bΓ=50.


Importantly, for constant parameters ω0, A, and B the solution of Eq. (14) can be found in textbooks (see, e.g., pages 86–87 in Landau and Lifshitz1982), and it has the form

(17) δ y = F ω 0 2 + a cos ψ - A a 2 2 ω 0 2 + A a 2 6 ω 0 2 cos ( 2 ψ ) + O ( a 3 ) ,

where ψ=ωt+ϕ0 (ϕ0 takes into account initial conditions), a is the amplitude of the linear part of δy, and

(18) ω = ω 0 + 3 B 8 ω 0 - 5 A 2 12 ω 0 3 a 2 .

Therefore, δy includes harmonic oscillations of amplitude a and O(a2) corrections and oscillation frequency ω (that includes an O(a2) correction to ω0). As is shown in the Appendix, when ω0 is a slow function of time as in our case [dω0/dtO(Γ)], the solution in Eq. (17) remains the same, but ψ is replaced by ψ=ωdt+ϕ0 and the oscillation's amplitude a becomes a slow function of time such that ωa2= const.


Figure 4Numerical solutions of x(t) in Eqs. (2)–(5) starting from the same initial conditions as in Fig. 3. The two black curves show the monotonic evolution of x described by Eq. (20) for t<tcr and by Eq. (21) for t>tcr, where tcr=50 as in Fig. 3. The black curves terminate near the Equator where the adiabaticity breaks down since the frequency, ω0, tends to 0 there according to Eq. (15).


This completes our solution for the latitude, y, and we proceed to the longitudinal dynamics. The dynamics in the zonal direction, x, are governed by Eq. (2), which after substitution of Eq. (13) becomes

(19) d x d t = D + y m 1 + b 2 y m + ( 1 + b y m ) δ y + b 2 δ y 2 .

Here again we consider two cases. For t<tcr, D+ym+(1+b2ym+)=0, and therefore by averaging locally in time over a single oscillation and using Eqs. (17) and (15), we get

(20) d x d t = - b a 2 2 + F ω 0 .

This equation shows that the average zonal drift is a nonlinear phenomenon in terms of the amplitude of oscillation. Since, as was shown above, F>0 for t<tcr, the drift is determined by the balance between ba2/2 and F/ω0=bΓ2/ω04. Thus, the sign (direction) of the zonal drift is independent of the sign of Γ.

For t>tcr, ym0=-1/b, so D+ym0(1+b2ym0)=Γt-1/(2b), and therefore

(21) d x d t = Γ t - 1 2 b + 1 4 b a 2 .

Figure 5The different water column trajectories for b=2 with initial conditions x=y=D=0 and V=0.002. (a) Large |Γ| (red curves), (b) small |Γ| (blue curves). The values of Γ are noted near each of the curves: thin curves denote negative (westward-directed) stresses, and thick curves denote positive (eastward-directed) stresses. The integration time is 2tcr in all cases, but the curves terminate just prior to reaching the Equator. Note that for Γ=0.00001 the integration time, 1bΓ, is 5000; i.e., the columns in the right panel complete several thousand oscillations in the course of the integration.


Figure 4 displays numerical solutions of x(t) in Eqs. (2)–(5) starting from the same initial conditions as in Fig. 3. The black curves show the monotonic evolution (averaged over oscillations) of x described by the theory developed here. Since the trajectories originate in midlatitudes, a westward-directed wind stress will always stir the trajectories away from the Equator, so Γ in Eq. (21) must be positive; i.e., the long-term zonal drift on the Equator has to be directed eastward.

Figure 5 compares the (x(t),y(t)) trajectories emanating from x=0=D=y and V=0.002 for two pairs of eastward-directed (thick curves) and westward-directed (thin curves) wind stresses. In the right panel the magnitude of the wind stress is small (0.0001), and in the left panel the magnitude of the wind stress is large (0.005). The two curves in each panel demonstrate that, as concluded above, the zonal drift is independent of the sign (direction) of the wind stress (since according to Eq. 20 it is proportional to F∝Γ2). A comparison between the trajectories in the two panels shows that for tiny wind stress (b) the trajectories drift westward as in the force-free, inertial oscillations, while with the increase in the magnitude of the wind stress (a) the zonal drift is directed eastward. In accordance with the intuits presented in the Introduction on the f-plane solution the oscillation's (inertial) frequency changes with latitude, i.e., increasing or decreasing in northward- or westward-directed trajectories, while the oscillation's amplitude follows the opposite pattern.

4 Discussion and summary

The two simple limits of b=0 (Ekman transport on the f plane) and Γ=0 (inertial trajectories on the β plane) should be discussed as special cases of the present theory. These limits are well known in physical oceanography, but they were never presented as limits of a single dynamical system.

In the b=0 limit (wind forced transport on the f plane) the potential in Eq. (8) becomes Φ(y)=12(D+y)2 (recall that Dt). This potential has a single minimum at yf=-D, and the frequency of oscillation near this point is ωf=1. Near yf the potential, Φ, is identical to that of the harmonic oscillator: 12(y-yf)2. Thus, the minimum yf=-D must decrease (or increase, depending on the sign of Γ) indefinitely at a rate that equals Γ; i.e., the potential Φ simply translates in the +y or y directions without changing its shape.

The Γ=0 limit (inertial trajectories on the β plane) implies, according to Eq. (4), that D is conserved, so the system in Eqs. (2)–(5) has two conserved quantities – D and the energy E. With the increase in the initial energy (say by increasing V(t=0)) the inertial trajectory will oscillate in (V,y) while drifting westward (see Ripa1997; Paldor2007) as on the sphere (Paldor2001). Equation (20) and the trajectories in Fig. 5 show that the long-term westward drift on the β plane when Γ≠0 is slower than in the inertial case, Γ=0, and it is independent of the sign of Γ.

The solutions of the nonlinear system in Eqs. (2)–(5) are determined by the two initial conditions V(t=0) and y(t=0) (recall that x(t=0)=0=D(t=0) can be assumed without loss of generality since x does not affect the dynamics and D can be translated in time) and the values of the two parameters, b and Γ (that represent the dimensional parameters β and τx, respectively), for a total of four parameters! Thus, these solutions display a range of temporal evolution, and this work describes and analyzes the general properties of these solutions and illustrates them in numerical examples. In particular, the zonal drift of the trajectories can be eastward (as in Fig. 1 and the left panel of Fig. 5) or westward (as in the right panel of Fig. 5). The sensitive dependence of the drift on parameter values (including initial conditions) is a defining property of nonlinear systems such as that studied here.

The intent of the analysis in this work is to provide an overview of the complex phenomena that result from the extension of Ekman's theory to the β plane. In particular, this work shows that the zonal drift is independent of the sign of τx but depends on a (previously unknown) balance between V(0)2 (or the initial displacement from (ym+)2) and (τx)2. The values of the parameters used in the numerical results presented here were chosen to highlight the phenomena being discussed while still being realistic. Thus, with the velocity scale of f0Re=640 m s−1 the value of V(0)=0.002 used in Figs. 35 corresponds to a dimensional velocity of about 1 m s−1. Trajectories of much higher oscillation amplitudes will be encountered with higher V(0) values.

The symmetry between ym+ and ym- in the present theory suggests that for the same wind stress, Γ, the Southern Hemisphere's fixed point will also move towards the Equator, i.e., northward. However, in all other respects the evolution near ym- is identical to that described above for ym+. The importance of latitudes at which the curl of the wind stress vanishes, which play a fundamental role in Stommel (1948) vorticity-based theory of wind-driven ocean gyres, cannot be captured in extensions of the present Lagrangian theory. However, extensions of the present new Lagrangian theory on the β plane can include variable zonal wind stress, τx(y), which can highlight the role played by latitudes at which the wind stress itself vanishes. The application of the concepts developed here to spherical geometry and to wind-driven circulation over the continental shelf (where the sloping bottom yields the topographic β effect) is an interesting goal for future studies.

Appendix A: Adiabatic evolution of meridional oscillations and initial conditions

In this Appendix we discuss adiabatic (slow) evolution of linear longitudinal oscillations described by (see Eq. 14)

(A1) d 2 δ y d t 2 = - ω 0 2 ( t ) δ y + F ω 0 2

and seek the solution of this equation of the form

(A2) δ y = a ( t ) cos 0 t ω 0 ( t ) d t + ϕ 0 + F ω 0 2 .

Here ϕ0 is added to take into account initial conditions, and we assume that the change in ω0 during one period 2π/ω0 of oscillations is small, i.e.,

(A3) d ω 0 d t 2 π ω 0 ω 0 ,

which is guaranteed if Γ is sufficiently small. This is our adiabaticity criterion. A similar condition, dadt2πω0a, is also assumed for the amplitude of oscillations. Next, we substitute Eq. (A2) into Eq. (A1) and neglect d2a/dt2 to get

(A4) 2 d a d t ω 0 + a d ω 0 d t = 0 ,


(A5) ω 0 a 2 = I = const .

The constant I (the action) is given by initial conditions. When the nonlinear terms in Eq. (14) are included in the analysis, the entire derivation of the weakly nonlinear solution as described in pages 86–87 of Landau and Lifshitz (1982) is not affected by the replacement of the linear component acos (ωt+ϕ0) by a(t)cos(0tω(t)dt+ϕ0) in the adiabatic problem, which is the basis of the solution in Eq. (17) in Sect. 3.

Finally, the action I, which remains constant all times, can be calculated from the initial conditions, δy(0) and V(0)=d(δy)/dt|t=0. Using Eq. (A2) we have δy(0)=a(0)cosϕ0+F(0)/ω02(0) and V(0)=-a(0)ω0(0)sinϕ0. Then

(A6) a ( 0 ) 2 = δ y ( 0 ) - F ( 0 ) ω 0 2 ( 0 ) 2 + V ( 0 ) ω 0 ( 0 ) 2


(A7) I = ω 0 ( 0 ) δ y ( 0 ) - F ( 0 ) ω 0 2 ( 0 ) 2 + V 2 ( 0 ) ω 0 ( 0 ) .

The case depicted in Figs. 3 and 4 has δy(0)=0, so one gets a2(0)=F2(0)+V2(0)ω02(0) and I=ω0(0)a2(0).

Data availability

No data sets were used in this article.

Author contributions

The research on the problem was initiated by NP, who also proposed the transformation to the pseudo-angular momentum, while LF proposed the application of the adiabaticity theory. Both authors contributed equally to the numerical calculations and paper preparation.

Competing interests

The contact author has declared that neither of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Review statement

This paper was edited by Anne Marie Tréguier and reviewed by Nicolas Grisouard and two anonymous referees.


Constantin, A. and Johnson, R.: Ekman-type solutions for shallow-water flows on a rotating sphere: a new perspective on a classical problem, Phys. Fluid., 31, 021401,, 2019. a

Ekman, V. W.: On the influence of the earth's rotation on ocean-currents., Ark. Mat. Astr. Fys., 2, 1–52, 1905. a

Gill, A. E.: Atmosphere-ocean dynamics, Vol. 30, Academic Press, ISBN 0122835204, 1982. a, b, c

Goldstein, H.: Classical mechanics, Eddison-Wesley Publishing Company, ISBN 978020131611, 1980. a

Landau, L. D. and Lifshitz, E. M.: Mechanics, Butterworth-Heinemann, Oxford, England, 3 Edn., ISBN 9780750628969, 1982. a, b

Paldor, N.: The zonal drift associated with time-dependent particle motion on the earth, Q. J. Roy. Meteorol. Soc., 127, 2435–2450,, 2001. a, b

Paldor, N.: The transport in the Ekman surface layer on the spherical Earth, J. Mar. Res., 60, 47–72,, 2002. a

Paldor, N.: Inertial particle dynamics on the rotating earth, in: Lagrangian Analysis and Prediction of Coastal and Ocean Dynamics, edited by: Griffa, A., Kirwan, A. D., Mariano, A. J., Ozgokmen, T., and Rossby, H. T., Cambridge University Press, ISBN 9780521870184, chap. 5, 119–135, 2007. a, b, c

Pedlosky, J.: Geophysical fluid dynamics, Vol. 710, Springer, ISBN 3540963871, 1987. a, b

Ripa, P.: “Inertial” Oscillations and the β-Plane Approximation (s), J. Phys. Oceanogr., 27, 633–647,<0633:IOATPA>2.0.CO;2, 1997. a, b

Stommel, H.: The westward intensification of wind-driven ocean currents, Eos, Trans. Am. Geophys. Union, 29, 202–206,, 1948. a

Vallis, G. K.: Atmospheric and oceanic fluid dynamics, Cambridge University Press, ISBN 9781107065505, 2017. a, b

Von Arx, W. S.: An introduction to physical oceanography, USA, ISBN 0201081741, 1964. a

Short summary
Mean surface transport forced by wind blowing over an ocean can assume many directions relative to the wind direction and does not have to be perpendicular to the direction of the wind. This is in contrast to a simpler 120-year-old theory that completely ignored Earth's sphericity and predicted that the mean ocean current should always be perpendicular to the direction of the overlying wind. In the new theory the direction of the mean current is determined by the values of several parameters.