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

. The seminal, Ekman (1905)’s, 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, 5 which 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 10 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.

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 25 wind-driven dynamics on the f -plane fully understood and quantified, the β-plane offers an in-between set-up where 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 suggest heuristically that the extension of Ekman's transport theory to the β-plane should include the following qualitative elements: 1. An increase/decrease in mean meridional velocity for an eastward/westward directed stress due to the decrease/increase 30 in Coriolis frequency when the water column moves southward/northward.
2. The frequency of oscillation about the mean velocity should decrease/increase (so oscillation period should increase/decrease) due to the decrease/increase in Coriolis frequency along the trajectory (for an eastward directed stress in the northern hemisphere while the opposite changes occur for westward directed tress and in the southern hemisphere).
3. Since the oscillation's frequency and amplitude are inversely correlated (energy flux is unchanged) a decrease in fre-35 quency 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). 40 The numerical solutions of the governing Lagrangian equations (see section 2 below) shown in figure 1 fully confirm the first 3 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 figure 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.

45
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 in a latitude where the Coriolis frequency vanishes -the equator. From that point onward the water column is subject to non-rotating dynamics 50 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 section 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 section 3 and the work concludes with a discussion and summary in section 4. 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 (Gill, 1982;Pedlosky, 1987;Vallis, 2017). 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 60 velocity components and the changes in the coordinate of the moving column i.e.: Here τ x is the uniform zonally directed wind stress (which is positive/negative for eastward/westward directed wind, respectively), ρ is the water density, H is the depth (thickness) of the layer, f = f 0 + βy is the Coriolis parameter (where f 0 = 2Ωsin(ϕ 0 ), β = 2Ωcos(ϕ 0 )/R e with R e and Ω -Earth's radius and rotation frequency, respectively and ϕ 0 -the latitude 65 where 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 details in e.g. chapter 9 of Gill (1982) is that here the Coriolis frequency, f , in the momentum equations is y-dependent.
The 4-dimensional system (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 1 f0 , x and y on R e so the velocity scale is f 0 R e . With this scaling the nondimensional Coriolis frequency is 1 + by where b = βRe f0 = cot(ϕ 0 ) is the nondimensional β. The system is further simplified by replacing U by the pseudo angular momentum, defined as D = U − y(1 + b 2 y) 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 75 a simple relation between the zonal velocity and the latitude (Paldor, 2001). 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 system (1) transforms to: Here t, x, y and V denote the nondimensional counterparts of the dimensional variables denoted by the same symbols in system (1) and, as explained above, D = U − y(1 + b 2 y) is the nondimensional pseudo angular momentum. Equation (4) (3) and (5) to the single second-order equation 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 equation (A3) in the Appendix]. We proceed by rewriting Eq. (6) as 105 Equation (7) describes the dynamics of a quasi-particle in a slowly (for small Γ) time varying quasi-potential well Φ(y, t). In figure 2 we illustrate this potential for Γ = b = 0.1 at times t = 0, 10, 20, ..., 100. The minima of these potentials, denoted collectively by y m , are given by the 3 roots of: Two cases should be considered depending on time being below or above the critical time For t < t cr , there exist two minima defined by Γt + y m 1 + 1 2 by m = 0, i.e., while for t > t cr , there exists a single minimum located at 115 Figure 3 shows the numerical solutions of y(t) when the column originates near y + m 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 y + m and y 0 m shown by the black curves. The direction of evolution of y + m (t) and its transformation to the constant y 0 m (t) at t = t cr correspond to the evolution of the red circles in figure 2 for t < t cr (left panel) and t > t cr (right panel). The averaged numerical solution is expected to deviate appreciably from the simple scenario shown here only for high 120 oscillation amplitude near y + m . 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, y + m , by the adiabatic theory (see pages 531-535 in Goldstein, 1980) it will stay near this minimum for t < t cr . At t = t cr , y + m transforms into y 0 m and therefore at all t > t cr , the system remains near y 0 m . Thus, the column remains near the minimum of Φ at all times, while this minimum is slowly decreasing for t < t cr and stays constant for t > t cr . Since the 125 trajectory originates near the minimum y + m 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 y = y m (t) + δy where y m (t) starts at y + m and later (i.e. at t = t cr ) transforms into y 0 m and δy is a small perturbation. We substitute this form of solution into Eq. (6) and rewrite the resulting equation as where F = −d 2 y m /dt 2 is an inhomogeneous forcing term and the coefficients of the other 3 terms on the RHS of this equation are:  (11) for t < tcr and ofy 0 m (t) given by equation (12) for t > tcr. The values of b = 2 and Γ = 0.001 used here imply that the change between the two approximate solutions occurs at tcr = 1 2bΓ = 50.
In the present model, equation (11) implies that for t < t cr , d 2 y + m /dt 2 = −bΓ 2 (1 − 2bΓt) −3/2 so according to equation (15) F = −d 2 y + m /dt 2 = bΓ 2 /ω 3 0 > 0 . The second term on the RHS of Eq. (14) describes linear oscillations having slowly varying 140 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 y m = y + m the term y + m (t) in Eq. (13) describes slow monotonic variation of the latitude shown by the black arrows in figure 2 at t < t cr . No such variation exists at t > t cr since then y m = y 0 m = const. As will be shown below, the nonlinear terms in (14) mostly affect the zonal drift in x.
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 145 Landau and Lifshitz, 1982) and it has the form where ψ = ωt + ϕ 0 , (ϕ 0 takes into account initial conditions), a is the amplitude of the linear part of the δy and Therefore, δy includes harmonic oscillations of amplitude a and O(a 2 ) corrections and oscillation frequency ω (that includes  (21) for t > tcr where tcr = 50 as in figure 3. The black curves terminate near the equator where the adiabaticity breaks down since the frequency, ω0, tends to 0 there according to equation (15) . This completes our solution for the latitude, y, and we proceed to the longitudinal dynamics. The dynamics in the zonal direction, x, is governed by Eq.
(2) which after substitution of (13) becomes Here again we consider two cases. For t < t cr , D + y + m (1 + b 2 y + m ) = 0 and, therefore, by averaging locally in time over a single oscillation and using Eqs. (17) and (15) we get This equation shows that the average zonal drift is a nonlinear phenomenon in terms of the amplitude of oscillation. Since, as 160 was shown above, F > 0 for t < t cr , the drift is determined by the balance between ba 2 /2 and F/ω 0 = bΓ 2 /ω 4 0 . Thus, the sign (direction) of the zonal drift is independent of the sign of Γ. 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, 1 bΓ , is 5000 i.e. the columns in the right panel complete several thousand oscillations in the course of the integration trajectories originate in mid-latitudes, a westward directed wind-stress will always stir the trajectories away from the equator so Γ in equation (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 170 is small (=0.0001) and in the left panel the magnitude of the wind tress 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 equation (20) it is proportional to F ∝ Γ 2 ). A comparison between the trajectories in the two panels shows that for tiny wind stress (right panel) the trajectories drift westward as in the force-free, inertial, oscillations while with the increase in the magnitude of the wind-stress (left panel) 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/decreasing in northward/westward directed trajectories while the oscillation's amplitude follow the opposite pattern.

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 180 presented as limits of a single dynamical system. In the b = 0 limit (wind forced transport on the f -plane) the potential in (8) becomes Φ(y) = 1 2 (D + y) 2 (recall: D = Γt). This potential has a single minimum at y f = −D and the frequency of oscillation near this point is ω f = 1. Near y f the potential, Φ, is identical to that of Harmonic Oscillator: 1 2 (y − y f ) 2 . Thus, the minimum y f = −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 185 without changing its shape.
The Γ = 0 limit (inertial trajectories on the β-plane) implies, according to equation (4), that D is conserved so system (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 Ripa, 1997;Paldor, 2007) as on the sphere (Paldor, 2001). Equation (20) and the trajectories shown in figure 5 show that the long-term westward drift on the β-plane when Γ ̸ = 0 190 is slower than in the inertial, Γ = 0, case and it is independent of the sign of Γ.
The solutions of the nonlinear system (2) -(5) are determined by the 2 initial conditions V (t = 0) and y(t = 0) (recall: 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 4 parameters! Thus, these solutions display a range of temporal evolution and this work describes and analyzes The intent of the analysis in this work is to provide an overview of the complex phenomena that result from the extension 200 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 (y + m ) 2 ) and (τ x ) 2 . The values of the parameters used in the numerical results presented here were chosen so as to highlight the phenomena being discussed while still being realistic. Thus, with the velocity scale of f 0 R e = 640m/s the value of V (0) = 0.002 used in figures 3 -5 corresponds to a dimensional velocity of about 1ms −1 . Trajectories of much higher oscillation amplitudes will be encountered 205 with higher V (0) values.
The symmetry between y + m and y − m 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 y − m is identical to that described above for y + m . The importance of latitudes where the curl of the wind-stress vanishes, that play a fundamental role in (Stommel, 1948) vorticity based theory of wind driven ocean gyres, can not be captured in extensions of the present 210 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 where 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.
Here ϕ 0 is added to take into account initial conditions and we assume that the change of ω 0 during one period 2π/ω 0 of 220 oscillations is small, i.e.
which is guaranteed if Γ is sufficiently small. This is our adiabaticity criterion. A similar condition, da dt 2π ω0 ≪ a, is also assumed for the amplitude of oscillations. Next, we substitute (A2) into (A1) and neglect d 2 a/dt 2 to get
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 manuscript preparation.