puntininonlineardynamicswinddrift2025.md 43 KB
category: literaturenote citekey: puntininonlineardynamicswinddrift2025 title: Nonlinear Dynamics of Wind-Drift Currents at Mid-Latitudes authors: "Puntini, Christian" year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.10454v1" zotero_key: PU9FBZ8M zotero_storage: 9P62JMZT collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Puntini, Christian"
status: converted
Nonlinear Dynamics of Wind-Drift Currents at Mid-Latitudes
Christian Puntini
Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090, Vienna, Austria
E-mail: christian.puntini@univie.ac.at
13 June 2025
Abstract. Starting from the Navier-Stokes equation in the f-plane approximation, we provide an exact and explicit solution of the governing equations at leading order for fluid flows in the upper layer of the ocean, driven by a wind stress. Such a solution highlights the presence of a mean Ekman current superimposed to trochoidal oscillations and a background geostrophic current.
Keywords: upper ocean; Ekman spiral; asymptotic analysis; Lagrangian description; Gerstner waves
Submitted to: Fluid Dyn. Res.
1. Introduction
The upper ocean, namely the topmost layer of sea water between the surface and the thermocline, with depths of the order of 50 m and rarely exceeding 100 m, is the region of the water column featuring the most complex dynamics. It is the most important stratum for the exchange of heat between the ocean and the atmosphere, as the first 2 − 3 m of the ocean posses the same heat capacity as the entire atmosphere above it. Moreover, around 50% of the solar radiation is absorbed within the top 0.5 m of the ocean. Other energetic processes occur in this upper layer, such as wave breaking: about half of the breaking surface wave kinetic energy is dissipated within 20% of the significant wave height from the surface (see Soloviev & Lukas 2014). Also, this water stratum is highly affected by meteorological events, such as rain, contributing to the surface and volume freshwater, momentum and heat fluxes. However, along with the physical processes cited above and many others (we refer to Soloviev & Lukas 2014 and Talley et al 2011 for a detailed description of the other phenomena occurring in the first meters of the ocean), there is another one particularly important: turbulence. In fact, the surface layer (also named surface mixing layer, as intensive mixing due to turbulence occurs there) is chaotic. Nevertheless, some coherent structures still emerge, such as wind-driven (Ekman) flows, and tidal and pressure-driven (geostrophic) currents. Ekman flows (or Ekman spirals), named after V. W. Ekman, who in 1905 was the first to provide a quantitative explanation for the phenomenon, describe the deflection of steady wind-driven ocean currents to the right of the prevailing wind direction in the Northern Hemisphere. This deflection arises from the momentum balance between the Coriolis acceleration, which acts on the upper ocean layers, and the frictional forces generated by turbulent stress induced by the wind (Ekman 1905). The deflection angle is to the left in the Southern Hemisphere.
Ekman's classical theory predicts a deflection angle of 45◦ ; however, in situ measurements have not only confirmed the existence of the Ekman spiral but also highlighted significant deviations between the observed deflection angles and those predicted by the classical solution (see Polton et al 2013, Roach et al 2015, and Yoshikawa & Masuda 2009).
In this paper, we extend the classical Ekman solution to investigate wind-driven nonlinear flows in the upper ocean layer, focusing on mid-latitude regions—areas away from the Equator and excluding the poles. A similar analysis for the Arctic Ocean has been conducted in Constantin 2022a, Constantin 2022b, and Puntini 2025. However, these solutions are not applicable to equatorial regions, where the Coriolis force vanishes at the Equator, and its sign reversal across the Equator generates azimuthal (i.e., longitudinal) flows (see, e.g., Constantin & Ivanov 2019 and Constantin & Johnson 2015).
Our approach involves the application of systematic theoretical methods to address problems in physical oceanography, following the framework established by works such as Constantin & Johnson 2016, Johnson 2022, and Johnson 2023. While it is undeniable that data are essential to support both new ideas and existing theories, we contend that valuable insights into the dynamics of geophysical flows can still be achieved through the careful (asymptotic) analysis of the equations of motion and the boundary conditions that define the problem. (see also Johnson 2017 and Johnson 2019 for a review of this methods), especially when studying near-inertial fluid flows, that seems to be hard to detect by satellite methods or via fixed moorins (see the discussion in Constantin & Monismith 2017).
The flows we are considering have length scales that are small compared to the dimensions of the Earth. In particular, as it will be detailed later, they are relevant for regions of about 1◦ of amplitude. Therefore, a reasonable first choice is to use local Cartesian coordinates and to write the Navier-Stokes equation in the f-plane approximation (Vallis 2017). However, we point out that the f-plane approximation has some limits, as is obtained by removing - in some sense - ad hoc terms from the general governing equations of fluid motion in spherical coordinates, and using tangent plane cartesian coordinates (see Constantin & Johnson 2019a, and Vallis 2017). The extension of this analysis to a more general setting will be done in future work.
We will present an exact solution to the nonlinear governing equations - at leading order - in Lagrangian framework, by specifying the trajectories of the individual fluid particles. These paths are composed of a mean Ekman spiral, trochoidal time-dependent near-inertial oscillations and a background geostrophic current.
For simplicity we will analyze the fluid motion in the Northern Hemisphere, but the results applies immediately in the Souther Hemisphere where the deflection angle is to the left.
2. Governing equations
We introduce a local Cartesian coordinate system with unit vectors (e1, e2, e3) and coordinates (x ′ , y′ , z′ ). The (x ′ , y′ , z′ )-system is associated with a point fixed on the sphere (other than at the two poles where e1 and e2 are not well-defined) which is rotating about its polar axis (with constant angular speed Ω′ ≈ 7.29 · 10−5 rad s−1 ). The associated velocity components are (u ′ , v′ , w′ ) =: u ′ . See Figure 1.
We use prime to denote physical (dimensional) quantities; the prime will be removed
Figure 1. The Cartesian coordinate system on the rotating Earth. Here, r represents the position vector of a point located at latitude θ and longitude φ. Z represents the null island.
after introducing the appropriate set of non-dimensional variables.
We adopt the f-plane approximation, as we are considering flows with small length scales (compared to Earth's dimensions), hence the Coriolis parameters f ′ = 2Ω′ sin θ and ˆf ′ = 2Ω′ cos θ can be considered constant, with θ denoting the fixed angle of latitude relative to the sector of the ocean being analyses. The local tangential coordinates x ′ and y ′ , pointing eastward and northward, respectively, will describe the position of points within the small regions under consideration, for every fixed z ′ . The coordinate z ′ describes the distance from the Earth's center in the radial direction. As the flows we are studying here are those present in the topmost stratum of the ocean, we anticipate that z ′ will be rewritten as z ′ = R′ + ζ ′ , where R′ ≈ 6371 km is the radius of the Earth (considered a sphere of constant geopotential).
The Navier-Stokes equation - in the f-plane approximation - is therefore written as (see Vallis 2017 or White 2002),
$$\frac{D\mathbf{u}'}{Dt'} + (\hat{f}'w' - f'v', f'u', -\hat{f}'u') =$$
$$= -\frac{1}{\rho'(z)}\nabla'p' + (0, 0, -g') + A'{h}\left(\frac{\partial^{2}}{\partial x'^{2}} + \frac{\partial^{2}}{\partial y'^{2}}\right)\mathbf{u}' + A'{v}\frac{\partial^{2}}{\partial z'^{2}}\mathbf{u}'$$ (1)
where A′ h and A′ v are the horizontal and vertical kinematic eddy viscosities, respectively; they are assumed to be constant, and we have taken the density to be only depthdependent, that is ρ ′ = ρ ′ (z), which is common for oceanic flows (see e.g. Constantin & Johnson 2023). In (1) we have used the common shorthand notation for the gradient and for the material derivative
$$\nabla' = \left(\frac{\partial}{\partial x'}, \frac{\partial}{\partial y'}, \frac{\partial}{\partial z'}\right), \quad \frac{D}{Dt'} = \frac{\partial}{\partial t'} + u'\frac{\partial}{\partial x'} + v'\frac{\partial}{\partial y'} + w'\frac{\partial}{\partial z'}, \tag{2}$$
where t ′ is the time variable. The equation of mass conservation is
$$\frac{D\rho'}{Dt'} + \rho' \left( \frac{\partial u'}{\partial x'} + \frac{\partial v'}{\partial y'} + \frac{\partial w'}{\partial z'} \right) = 0, \tag{3}$$
which, due to the assumption on the viscosity (implying Dρ′ Dt′ = w ′ dρ′ dz′) reduces to
$$\frac{w'}{\rho'}\frac{d\rho'}{dz'} + \frac{\partial u'}{\partial x'} + \frac{\partial v'}{\partial y'} + \frac{\partial w'}{\partial z'} = 0. \tag{4}$$
As anticipated, we write z ′ = R′ + ζ ′ , with R′ being the Earth's radius. We simplify the Navier-Stokes equation (1) by redefinig the pressure as
$$p' = P{\text{atm}} + g' \int{z'}^{R'} \rho'(r) dr + P'(x', y', z', t'),$$ (5)
giving
$$\frac{D\boldsymbol{u}'}{Dt'} + (\hat{f}'\boldsymbol{w}' - f'\boldsymbol{v}', f'\boldsymbol{u}', -\hat{f}'\boldsymbol{u}') = -\frac{1}{\rho'(z)}\nabla'P' + A'{h}\left(\frac{\partial^{2}}{\partial x'^{2}} + \frac{\partial^{2}}{\partial y'^{2}}\right)\boldsymbol{u}' + A'{v}\frac{\partial^{2}}{\partial z'^{2}}\boldsymbol{u}'. \tag{6}$$
Now we non-dimensionalise the equations of motion (6) and (4) by introducing new variables x, y, z, u, v, w, ρ, P according to the transformation
$$(x',y') = L'(x,y), \quad \zeta' = D'z, \quad \rho' = \overline{\rho'}\rho, \quad P' = \overline{\rho'}U'^2P,$$
$$(u',v',w') = U'(u,v,kw), \quad t' = \frac{L'}{U'}t.$$
$$(7)$$
Here L ′ is a length scale, D′ is a suitable depth scale (to be chosen later), U ′ an appropriate speed scale, ρ ′ is the average density and k a scaling factor.
We assume the length scale L ′ to be of the order - but not in excess - of the Rossby radius of deformation L ′ (θ). We report here the Rossby radius at four latitudes (see Chelton et al 1998)
$$\mathcal{L}'(10^{\circ}) \approx 100 \,\text{km}, \quad \mathcal{L}'(20^{\circ}) \approx 50 \,\text{km}, \quad \mathcal{L}'(45^{\circ}) \approx 30 \,\text{km}, \quad \mathcal{L}'(60^{\circ}) \approx 20 \,\text{km}.$$ (8)
We introduce the "thin-shell" parameter
$$\varepsilon = \frac{D'}{R'},\tag{9}$$
and we denote by K the ratio between L ′ and R′ ,
$$\mathfrak{K} = \frac{L'}{R'}. (10)$$
The asymptotic structure will be driven by the parameter ε. Because the length scales we are considering are small (see (8)), the regions we are considering will not exceed 1 ◦ of latitude, as 100 km ≈ 53 nmi and 60 nmi = 1◦ of latitude. Here nmi stands for "nautical miles". As a consequence, the f-plane approximation emerges as a suitable initial choice for the study of such problems, even if it possesses some limitations (see the discussion in Constantin & Johnson 2019a). We hope, in a subsequent work, to extend this kind of analysis to a more general setting.
After the nondimensionalisation (7), the Navier-Stokes equation (6) reads as (in component form)
$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + k \frac{L'}{D'} w \frac{\partial u}{\partial z} - fv + fkw \cot \theta = = -\frac{1}{\rho(z)} \frac{\partial P}{\partial x} + \frac{L'^2}{D'^2} \frac{1}{\text{Re}_v} \frac{\partial^2 u}{\partial z^2} + \frac{1}{\text{Re}_h} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} - \tan \theta \frac{\partial u}{\partial y} \right), \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + k \frac{L'}{D'} w \frac{\partial v}{\partial z} + fu = = -\frac{1}{\rho(z)} \frac{\partial P}{\partial y} + \frac{L'^2}{D'^2} \frac{1}{\text{Re}_v} \frac{\partial^2 v}{\partial z^2} + \frac{1}{\text{Re}_h} \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} - \tan \theta \frac{\partial v}{\partial y} \right), k \left( \frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + k \frac{L'}{D'} w \frac{\partial w}{\partial z} \right) - fu \cot \theta = = -\frac{1}{\rho(z)} \frac{L'}{D'} \frac{\partial P}{\partial z} + \frac{L'^2}{D'^2} \frac{k}{\text{Re}_v} \frac{\partial^2 w}{\partial z^2} + \frac{k}{\text{Re}_h} \left( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} - \tan \theta \frac{\partial w}{\partial y} \right),$$ (11)
where we have introduced the parameters
$$f = 2\frac{\Omega'L'}{U'}\sin\theta$$ , $\operatorname{Re}{h} = \frac{U'L'}{A'{h}}$ , $\operatorname{Re}{v} = \frac{U'L'}{A'{v}}$ . (12)
The first of these is an inverse Rossby number and the other two are Reynolds numbers. We have that f ≈ 25, if we use L ′ = L ′ (θ) at the specific latitudes (8).
The equation of mass conservation becomes
$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + k \frac{L'}{D'} \left( \frac{w}{\rho} \frac{d\rho}{dz} + \frac{\partial w}{\partial z} \right) = 0.$$ (13)
Following Constantin & Johnson 2019a, we make use of L′2 D′2Rev for defining the relevant depth-scale D′ by setting:
$$\frac{L'^2}{D'^2} = \text{Re}{\text{v}} \Longrightarrow D' = \sqrt{\frac{L'A'{\text{v}}}{U'}},\tag{14}$$
giving
$$\frac{1}{\mathrm{Re_h}} = \frac{D'^2}{L'^2} \mathfrak{X} = \frac{\varepsilon^2}{\mathfrak{K}^2} \mathfrak{X},\tag{15}$$
where
$$\mathfrak{X} := \frac{A{\rm h}'}{A{\rm v}'} = \frac{\mathrm{Re}{\rm v}}{\mathrm{Re}{\rm h}}.\tag{16}$$
For reference, we set the value of the mean vertical kinetic eddy viscosity (also referred to as the coefficient of vertical turbulent exchange, see Kamenkovich 1977) to be
$$A_{\rm v}' \approx 0.01 \,\rm m^2 s^{-1}$$ . (17)
Assuming U ′ ≈ 0.1 ms−1 , which is a typical value of ocean currents velocities (see e.g. Constantin & Johnson 2019b), and considering (8) and (17), we find that the values of D′ are in the range
$$45 \,\mathrm{m} \lesssim D' \lesssim 100 \,\mathrm{m},\tag{18}$$
which are, in fact, the standard values of the near-surface layer of the ocean, where wind-drift effects are significant (see Kamenkovich 1977), validating the choice of the depth scale D′ in (14).
We elect to consider the class of problems for which, in the limit ε → 0,
$$k = \varepsilon \kappa, \tag{19}$$
with the choice of κ = o(1), relevant for the problem we aim to study (see Constantin & Johnson 2019a). As an example, a possible definition for κ would be κ = ε n , with n > 0, as done in Constantin & Johnson 2024 in the context of atmospheric flows. However, as we are interested only in the leading order terms of the governing equations, the exact definition of κ is not important; it is sufficient to say κ = o(1).
Equations (11) and (13) therefore becomes
$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + \kappa \Re w \frac{\partial u}{\partial z} - fv + f\kappa \varepsilon w \cot \theta = = -\frac{1}{\rho(z)} \frac{\partial P}{\partial x} + \frac{\partial^2 u}{\partial z^2} + \frac{\varepsilon^2}{\Re^2} \Re \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} - \tan \theta \frac{\partial u}{\partial y} \right), \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + \kappa \Re w \frac{\partial v}{\partial z} + fu = = -\frac{1}{\rho(z)} \frac{\partial P}{\partial y} + \frac{\partial^2 v}{\partial z^2} + \frac{\varepsilon^2}{\Re^2} \Re \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} - \tan \theta \frac{\partial v}{\partial y} \right), \frac{\kappa \varepsilon^2}{\Re} \left( \frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + \kappa \Re w \frac{\partial w}{\partial z} \right) - \frac{\varepsilon}{\Re} fu \cot \theta = = -\frac{1}{\rho(z)} \frac{\partial P}{\partial z} + \frac{\kappa \varepsilon^2}{\Re^2} \frac{\partial^2 w}{\partial z^2} + \frac{\kappa \varepsilon^4}{\Re^3} \Re \left( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} - \tan \theta \frac{\partial w}{\partial y} \right),$$
and
$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \kappa \mathfrak{K} \left( \frac{w}{\rho} \frac{d\rho}{dz} + \frac{\partial w}{\partial z} \right) = 0. \tag{21}$$
At leading order - that is, taking the limit ε → 0, with κ = o(1) (see Constantin & Johnson 2019a) and keeping the other parameters fixed - the momentum equations (20) read as
$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} - fv = -\frac{\partial P}{\partial x} + \frac{\partial^2 u}{\partial z^2}, \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + fu = -\frac{\partial P}{\partial y} + \frac{\partial^2 v}{\partial z^2}, \frac{\partial p}{\partial z} = 0.$$ (22)
We remark the fact that the limit ε → 0 is not a real limit in the mathematical sense, but rather a formal way to obtain the equations at leading order.
The continuity equation (21) at leading order simplifies as
$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0. {23}$$
From (22) and (23), the flow appear to be a depth-dependent 2D horizontal flow.
2.1. Boundary conditions
For the full characterisation of the geophysical flow under consideration, a major role is played by the boundary conditions. We impose the following boundary conditions at the free surface ζ ′ = h ′ (x ′ , y′ , t′ ) (recall that we are writing z ′ = R′ + ζ ′ = R′ + D′ z):
• dynamic boundary condition: we impose a surface pressure, namely
$$p' = P_{\text{atm}} \qquad on \qquad \zeta' = h';$$ (24)
• kinematic boundary condition: we assume that no fluid particles will leave the fluid via the free surface, namely
$$w' = \frac{\partial h'}{\partial t'} + u' \frac{\partial h'}{\partial x'} + v' \frac{\partial h'}{\partial y'} \qquad on \qquad \zeta' = h'; \tag{25}$$
• shear stress: we require the shear stress at the surface of the water to be given. We express this in the form
$$\tau_1'(x', y', t') = \rho'(0)A_v'\frac{\partial u'}{\partial z'}, \quad \tau_2'(x', y', t') = \rho'(0)A_v'\frac{\partial v'}{\partial z'}, \qquad on \qquad \zeta' = h'$$ (26)
Moreover, we require that the water is in hydrostatic state at great depth. This boundary condition will be detailed later.
The boundary conditions also need to be non-dimensionalised according to (7), and, analogously to the equations of motion, studied in the limit ε → 0 to retain the leading order terms.
Starting from the kinematic boundary condition (25)
$$\kappa \varepsilon w = \frac{H'}{L'} \left( \frac{\partial h}{\partial t} + u \frac{\partial h}{\partial x} + v \frac{\partial h}{\partial y} \right) \qquad on \qquad \zeta' = h', \tag{27}$$
where we have written h ′ = H′h. As the depth scale H′ is independent of ε, in the limit ε → 0, (27) reads as
$$\frac{\partial h}{\partial t} + u \frac{\partial h}{\partial x} + v \frac{\partial h}{\partial y} = 0 \qquad on \qquad z = \frac{H'}{D'} h. \tag{28}$$
The condition (28) implies the free surface is invariant for the flow, so the most reasonable physical assumption is that the free surface is flat. We can assume, without loss of generality, that the flat free surface is given by z = 0. Note that the presence of a flat free surface - namely the absence of surface waves - emerges as a consequence of the asymptotic behavior of the governing equations due to the scaling (7) relevant to the problem we are considering.
In particular, surface waves arise if u ′ and w ′ are of the same order (see Soloviev & Lukas 2014), while a fundamental assumption of the flow we are considering is that the vertical velocity w ′ has order of magnitude much lower that the horizontal velocities u ′ and v ′ , by an order of κε.
Due to the approximation of the free surface as a flat one, the condition (24) became
$$P = 0 \quad on \quad z = 0, \tag{29}$$
while, as the vertical velocity component is neglected at leading order in (22) - (23), the boundary conditions capturing the decay with depth of wind drift flows is given by
$$\frac{\partial (u'^2 + v'^2)}{\partial z'} < 0 \quad for \quad z' < 0, \tag{30}$$
which become
$$\frac{\partial (u^2 + v^2)}{\partial z} < 0 \quad for \quad z < 0, \tag{31}$$
when adimensionalised.
The stress generated by the wind blowing over the water is given by the bulk formula
$$\tau_a' = \rhoa' C{aw} | \boldsymbol{u}_a' - \boldsymbol{u}_s' | (\boldsymbol{u}_a' - \boldsymbol{u}_s'), \tag{32}$$
where ρ ′ a is the air density, Caw is the (nondimensional) drag coefficient of air over water and u ′ a = (u ′ a , v′ a ) is the wind velocity describing the horizontal momentum exchange at the surface and u ′ s = (u ′ s , v′ s ) is the velocity field of the ocean at the surface. Namely, if there is a current present at the sea surface, as one could expect, the only contribution to the stress generated by the wind is given by the velocity of the wind relative to the ocean surface current, and not just by the wind velocity in its entirety. The wind velocity is usually referred as the wind velocity at 10 m above the sea surface. For order of magnitude estimates, we follow Soloviev & Lukas 2014, and set
$$C_{aw} \approx \begin{cases} 1.1 \cdot 10^{-3}, & for \quad |\mathbf{u}'_a| \lesssim 5 \,\mathrm{m \, s^{-1}}, \ (0.61 + 0.063 |\mathbf{u}'_a|) \cdot 10^{-3}, & for \quad 5 \,\mathrm{m \, s^{-1}} \lesssim |\mathbf{u}'_a| \lesssim 22 \,\mathrm{m \, s^{-1}}, \end{cases}$$ (33)
the second one being the linear relationship proposed by Smith 1988. For winds blowing over 22 m s−1 , a formula for the drag coefficient is difficult to obtain (see Soloviev & Lukas 2014). These are extreme wind conditions that eventually occur rarely and for short periods (usually these are gusts). As we are considering a time scale of the order of some days (see (7) and (8)), wind gusts over short time-scales are not relevant for our analysis. See Figure 2 for order of magnitude of the blowing wind at different latitudes.
Equaling (26) and (32), namely asking the continuity of the shear stress across the
Figure 2. Monthly-averaged ocean wind speed and direction vectors, with vector lengths proportional to the reference scale (in m s−1 ), based on observations from NASA's QuikSCAT satellite. Image credit: NOAA.
water surface, gives
$$\rho'(0)A'{v}\left(\frac{\partial u'}{\partial z'}, \frac{\partial v'}{\partial z'}\right) = \rho'{a}C{aw}\sqrt{(u'{a} - u'{s})^{2} + (v'{a} - v'{s})^{2}}\left(u'{a} - u'{s}, v'{a} - v'_{s}\right) \qquad on \qquad \zeta' = h'.(34)$$
Writing (u ′ a , v′ a ) = U ′ (ua, va) and (u ′ s , v′ s ) = U ′ (us, vs), the nondimensional form of (34) reads as
$$\left(\frac{\partial u}{\partial z}, \frac{\partial v}{\partial z}\right) = \gamma{aw} C{aw} \sqrt{(u_a - u_s)^2 + (v_a - v_s)^2} \left(u_a - u_s, v_a - v_s\right) \qquad on \qquad z = 0, (35)$$
where we have defined the adimensional parameter
$$\gamma_{aw} = \frac{\rho_a' U' D'}{\rho'(0) A_v'}.$$ (36)
3. The nonlinear solution in Lagrangian framework
The presence of a pressure gradient is offset by a geostrophic current, consistent with the ocean being in geostrophic balance (see Kundu 1990 and Roberti 2025).
Therefore, we decompose the horizontal velocity field in its geostrophic (pedices G) and ageostrophic components (b symbol) as
$$u = U_G + \widehat{u}, \qquad v = V_G + \widehat{v},$$ (37)
where the geostrophic velocity components satisfy
$$\frac{\partial P}{\partial x} = f V_G, \qquad \frac{\partial P}{\partial y} = -f U_G.$$ (38)
We assume UG and VG to be constant. Due to (37) and (38), the first two equations of (22), reduces to
$$\frac{Du}{Dt} - f\widehat{v} = \frac{\partial^2 u}{\partial z^2}, \qquad \frac{Dv}{Dt} + f\widehat{u} = \frac{\partial^2 v}{\partial z^2}, \tag{39}$$
while the third one, namely ∂P ∂z = 0 remains unchanged.
Adopting the Lagrangian approach (that is describing the motion of every single particle of fluid), and specifying at every time t the positions
$$x(t; a, b, z) = a + [d_1(z) + U_G] t - \frac{1}{k} e^{k(b+z)} \sin(k(a-z-ct)),$$
$$y(t; a, b, z) = b + [d_2(z) + V_G] t + \frac{1}{k} e^{k(b+z)} \cos(k(a-z-ct)),$$ (40)
of the horizontally moving fluid particles in terms of the depth z, the material variables a, b and the parameters k > 0 and c > 0, for suitably chosen functions d1(z) and d2(z), we will show that a solution to the nonlinear ice-drift flow at the North Pole (22), (23) satisfying the boundary conditions is given by (40). For the labeling variables, we have that a ∈ R, while b ∈ (b1, b0), with b1 < b0 < 0, ensuring an exponential decay for z ≤ 0 in (40). Writing
$$\xi = k(a - z - ct),\tag{41}$$
it follows that at every fixed time t we have
$$\begin{pmatrix} \frac{\partial x}{\partial a} & \frac{\partial x}{\partial b} \ \frac{\partial y}{\partial a} & \frac{\partial y}{\partial b} \end{pmatrix} = \begin{pmatrix} 1 - e^{k(b+z)}\cos(\xi) & -e^{k(b+z)}\sin(\xi) \ -e^{k(b+z)}\sin(\xi) & 1 + e^{k(b+z)}\cos(\xi) \end{pmatrix},\tag{42}$$
with inverse given by
$$\begin{pmatrix} \frac{\partial a}{\partial x} & \frac{\partial a}{\partial y} \ \frac{\partial b}{\partial x} & \frac{\partial b}{\partial y} \end{pmatrix} = \frac{1}{1 - e^{2k(b+z)}} \begin{pmatrix} 1 + e^{k(b+z)}\cos(\xi) & e^{k(b+z)}\sin(\xi) \ e^{k(b+z)}\sin(\xi) & 1 - e^{k(b+z)}\cos(\xi) \end{pmatrix}. \tag{43}$$
The determinant of the matrix in (42) is equal to 1 − e 2k(b+z) , hence the flow is areapreserving, implying that the continuity equation (23) is satisfied (see Constantin 2011). From (40) we have that
$$u = \frac{\partial x}{\partial t} = d_1(z) + U_G + ce^{k(b+z)}\cos(\xi), \tag{44}$$
$$v = \frac{\partial y}{\partial t} = d_2(z) + V_G + ce^{k(b+z)}\sin(\xi), \tag{45}$$
giving, for the ageostrophic components of the velocity field,
$$\widehat{u} = d_1(z) + ce^{k(b+z)}\cos(\xi),\tag{46}$$
$$\widehat{v} = d_2(z) + ce^{k(b+z)}\sin(\xi). \tag{47}$$
Moreover, since
$$\begin{pmatrix} \frac{\partial u}{\partial x} & \frac{\partial v}{\partial x} \ \frac{\partial u}{\partial y} & \frac{\partial v}{\partial y} \end{pmatrix} = \begin{pmatrix} \frac{\partial a}{\partial x} & \frac{\partial b}{\partial x} \ \frac{\partial a}{\partial y} & \frac{\partial b}{\partial y} \end{pmatrix} \begin{pmatrix} \frac{\partial u}{\partial a} & \frac{\partial v}{\partial a} \ \frac{\partial u}{\partial a} & \frac{\partial v}{\partial b} \end{pmatrix} \ = \frac{kc}{1 - e^{2k(b+z)}} \begin{pmatrix} -e^{k(b+z)}\sin(\xi) & e^{k(b+z)}\cos(\xi) + e^{2k(b+z)} \\ne^{k(b+z)}\cos(\xi) - e^{2k(b+z)} & e^{k(b+z)}\sin(\xi) \end{pmatrix}, (48)$$
at each horizontal level the two-dimensional horizontal flow (40) has constant vorticity
$$\frac{\partial u}{\partial y} - \frac{\partial v}{\partial x} = \frac{2kc}{1 - e^{-2k(b+z)}}. (49)$$
Using the expressions
$$\frac{Du}{Dt} = \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = kc^2 e^{k(b+z)} \sin(\xi), \tag{50}$$
$$\frac{Dv}{Dt} = \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = -kc^2 e^{k(b+z)}\cos(\xi), \tag{51}$$
from (39) we get the system of differential equations
$$fd_2(z) + d_1''(z) + (fc + 2k^2c - kc^2)e^{k(b+z)}\sin(\xi) = 0,$$ (52)
$$fd_1(z) - d_2''(z) + (fc + 2k^2c - kc^2)e^{k(b+z)}\cos(\xi) = 0.$$ (53)
Given the dispersion relation
$$c = \frac{f}{k} + 2k,\tag{54}$$
and writing in complex-variable notation
$$d(z) = d_1(z) + id_2(z), (55)$$
the system (52) reduces to
$$d''(z) - ifd(z) = 0, (56)$$
whose general solution is given by
$$d(z) = A e^{(1+i)\lambda z} + B e^{-(1+i)\lambda z},$$ (57)
where λ = qf 2 and A, B ∈ C.
As wind drift currents are insignificant at great depth, namely imposing the boundary condition (31), we have that B = 0. The solution of (56) therefore takes the form
$$d(z) = d(0)e^{(1+i)\lambda z}, (58)$$
with d(0) = d1(0) + id2(0) ∈ C. The time-average of the velocity (44) over a period T = 2π ck = 2π f+2k , which in view of (58) reads as
$$u + iv = d(0)e^{(1+i)\lambda z} + ce^{k(b+z)}e^{ik(a-z-ct)} + \mathbf{U}_G,$$ (59)
where UG = UG + iVG, yields the mean drift current
$$\langle u + iv \rangle_T = d(0)e^{(1+i)\lambda z} + \mathbf{U}_G. \tag{60}$$
Here ⟨·⟩T represents the time-average, and from (60) it is evident that d(0) is the mean wind-drift ocean current at the surface.
Due to (60), the boundary condition (35) reduces - in an averaged sense - to
$$\frac{\lambda}{C{aw}\gamma{aw}}(\mathfrak{d}_1 - \mathfrak{d}_2, \mathfrak{d}_1 + \mathfrak{d}_2) = |(u_a - u_s, v_a - v_s)|(u_a - u_s, v_a - v_s), \tag{61}$$
in which the unknowns are d1 = d1(0) and d2 = d2(0). Given that d0 is the mean wind-drift ocean current at the surface, we set us = UG + d1(0), vs = VG + d2(0), and (61) can be rewritten as (in complex notation and defining β = Caw√ γaw 2λ )
$$(\mathfrak{d}_1+i\mathfrak{d}_2)=e^{-i\frac{\pi}{4}}\beta\left|\left(u_a-U_G-\mathfrak{d}_1\right)+i(v_a-V_G-\mathfrak{d}_2)\right|\left[\left(u_a-U_G-\mathfrak{d}_1\right)+i(v_a-V_G-\mathfrak{d}_2)\right].(62)$$
Writing
$$(u_a - U_G - \mathfrak{d}_1) + i(v_a - V_G - \mathfrak{d}_2) = \mathcal{U}e^{i\Phi}, \tag{63}$$
the boundary condition (62) reads as
$$\mathcal{U}^2 e^{i\Phi} + \frac{1+i}{\beta} \mathcal{U} e^{i\Phi} = \frac{Z}{\beta},\tag{64}$$
where
$$Z = (u_a - U_G - v_a + V_G) + i(u_a - U_G + v_a - V_G).$$ (65)
For Z ̸= 0, we get from (64),
$$e^{i\Phi} = \frac{Z}{\mathcal{U}(\beta\mathcal{U} + (1+i))}. (66)$$
Multypling (66) by its complex conjugates gives the polynomial expression
$$\mathcal{P}(\mathcal{U}) = \beta^2 \mathcal{U}^4 + 2\beta \mathcal{U}^3 + 2\mathcal{U}^2 - |Z|^2 = 0.$$ (67)
As P(0) < 0, limR→±∞ P(R) = +∞ and P ′′(R) > 0, P is strictly convex, ensuring that there is one solution U > 0 of (67). Therefore, we point out that such computations have highlighted how the boundary condition (35) determines in a unique way the winddrift surface current d(0), knowing the geostrophic current and the wind speed. More precisely, due to (63) and (65), it is necessary just to know ua − UG and va − VG to determine the wind-drift current at the surface, d(0).
Writing the speed of the wind relative to the geostrophic current as
$$(u_a - U_G) + i(v_a - V_G) = We^{i\Theta}$$ (68)
and
$$d(0) = \mathfrak{d}_1 + i\mathfrak{d}_2 = \mathcal{D}e^{i\Psi} \tag{69}$$
we can rewrite the expression (62) as
$$(i+i)\mathcal{D}e^{i\Psi} = \sqrt{2}\beta \left| \mathcal{W}e^{i\Theta} - \mathcal{D}e^{i\Psi} \right| \left( \mathcal{W}e^{i\Theta} - \mathcal{D}e^{i\Psi} \right)$$ (70)
giving
$$\frac{(i+i)\mathcal{D}}{\sqrt{2}\beta \left| \mathcal{W}e^{i\Theta} - \mathcal{D}e^{i\Psi} \right|} + \mathcal{D} = \mathcal{W}e^{i(\Theta - \Psi)}. \tag{71}$$
The angle Θ − Ψ represents the deflection angle between the wind velocity (relative to the oceanic geostrophic current) and the wind-drift flow at the surface. As the left side of (71) is a complex number with imaginary part positive and smaller than the real part, it follows that 0 < Θ − Ψ < π/4. See Figure 3
Figure 3. Schematic depiction of the surface currents: in red the wind velocity; in blue the geostrophic current, in black the surface current with the angle Θ − Ψ to the right of the vector of the wind velocity relative to the geostrophic current (olive green)
Of particular interest is also the structure of the solution (40). In fact, for any given k > 0, the horizontal particle path given by (40) coupled with (58)
$$x(t;a,b,z) + iy(t;a,b,z) = (a+ib) + \left[d(0)e^{(1+i)\lambda z} + \mathbf{U}_G\right]t + \frac{1}{k}e^{k(b+z)}e^{i\left(\frac{\pi}{2} + k(a-z) - (f+2k^2)t\right)},$$ (72)
is a solution of (22) and (23) satisfying the boundary conditions and representing parametric trochoids. The mean ageostrophic current, over a period T = 2π ck = 2π f+2k 2 ,
$$\langle \widehat{u} + i\widehat{v} \rangle_T = \frac{1}{T} \int_0^T \left( \widehat{u}(t; a, b, z) + i\widehat{v}(t; a, b, z) \right) dt = d(0)e^{(1+i)\lambda z}$$ (73)
is the classical Ekman spiral (see Figure 4), while the oscillatory perturbation of the current
$$\frac{1}{k}e^{k(b+z)}e^{i\left(\frac{\pi}{2}+k(a-z)-(f+2k^2)t\right)}\tag{74}$$
is near-inertial, as f ≈ 15 and k << 1 gives
$$f + 2k^2 \approx f \tag{75}$$
for the frequency, hence the period of oscillations is approximately the inertial period of the Earth.
Figure 4. The wind-drift current. In blue the geostrophic current, in black the winddriven surface current, and in red the Ekman spiral. The u, v-axes, and the z-axis are not to scale.
Averaging over the depth the expression in (73) gives the depth-averaged mean-drift Ekman current of the solution (72), also known as (time-averaged) Ekman trasport,
$$\mathcal{I}{Ek} := \int{-\infty}^{0} d(0) e^{(1+i)\lambda z} dz = \frac{d(0)}{\lambda \sqrt{2}} e^{-i\frac{\pi}{4}}, \tag{76}$$
which is at 45◦ to the right of the mean surface current d(0); this is a classical result of Ekman's theory.
4. Discussion
The presented Lagrangian analysis provides an explicit solution for the nonlinear governing equations - at leading order - for the wind-driven flows in regions away from the poles (where classical spherical coordinates fails, and the use of rotated spherical coordinates is necessary - see Constantin & Johnson 2023 and Puntini 2025) and the equatorial band (due to the breakdown of the Ekman flow - see Constantin & Johnson 2015). The solution is the sum of a mean wind-drift current (the well-known Ekman spiral), a background geostrophic current balancing the pressure gradient, and nearinertial oscillations. For every fixed value of z, the particle path described (72) represents a trochoid, the only known solution for the 2D incompressible Euler equations with free boundary discovered by Gerstner (see e.g. Milne-Thomson 1938, Henry 2008, Constantin & Strauss 2010 or Abrashkin and Pelinovsky 2022). This study extends the solutions of Constantin 2022a, 2022b and Puntini 2025 for the Arctic Ocean, where ice and the Transpolar Drift Current are present, to mid-latitudes regions.
Recalling the momentum equations at leading order
$$\frac{Du}{Dt} - fv = -\frac{\partial P}{\partial x} + \frac{\partial^2 u}{\partial z^2},\tag{77}$$
$$\frac{Dv}{Dt} + fu = -\frac{\partial P}{\partial y} + \frac{\partial^2 v}{\partial z^2},\tag{78}$$
we can think the solution (72) as a nonlinear "superposition" of two terms, namely the Ekman spiral, solution of
$$-fv = \frac{\partial^2 u}{\partial z^2},\tag{79}$$
$$fu = \frac{\partial^2 v}{\partial z^2},\tag{80}$$
and the Gerstner's trocoidal solution of the Euler equations
$$\frac{Du}{Dt} = -\frac{\partial P}{\partial x},\tag{81}$$
$$\frac{Dv}{Dt} = -\frac{\partial P}{\partial y}.\tag{82}$$
We also have shown that, knowing the velocity of the wind and of the background geostrophic current, the wind-drift surface current is uniquely determined by the boundary conditions. Moreover, the deflection angle between the wind-drift flow at the surface and the blowing wind (relatively to the geostrophic flow) is less than 45◦ . In this sense, we also point out a the importance of the boundary conditions. In some textbooks or papers (see e.g. Trenberth 1990 or Marshall & Plumb 2016) the stress of the wind over the ocean is expresses as
$$\tau_a' = \rhoa' C{aw} | \boldsymbol{u}_a' | \boldsymbol{u}_a', \tag{83}$$
without taking into account the surface motion of the ocean. Assuming (83), equation (35) would be
$$\left(\frac{\partial u}{\partial z}, \frac{\partial v}{\partial z}\right) = \gamma{aw} C{aw} |(u_a, v_a)| (u_a, v_a), \tag{84}$$
which, due to (60), reads as
$$\frac{\lambda}{C{aw}\gamma{aw}}\left(\mathfrak{d}{1}-\mathfrak{d}{2},\mathfrak{d}{1}+\mathfrak{d}{2}\right)=\left|\left(u{a},v{a}\right)\right|\left(u{a},v{a}\right),\tag{85}$$
or equivalently as
$$(\mathfrak{d}_1 + i\mathfrak{d}_2) = e^{-i\frac{\pi}{4}}\beta\sqrt{u_a^2 + v_a^2}(u_a + iv_a),\tag{86}$$
where we have written again d1 = d1(0), d2 = d2(0) and β = Caw√ γaw 2λ . In this case, the deflection angle would be of exactly 45◦
We have assumed the vertical eddy viscosity to be constant, which is, to a first approximation, a good choice for oceanic flows; however, a direction that is of interest is to extend this model to the case of depth-dependent or time-dependent eddy viscosities, expanding the recent investigations by Bressan & Constantin 2019, Constantin et al 2020, Dritschel et al 2020, Constantin 2021, Roberti 2021 and Roberti 2022).
Acknowledgments
- The author declares no conflict of interest.
- The author is supported by the Austrian Science Fund (FWF) [grant number Z 387-N].
References
- Abrashkin A A, and Pelinovsky E N Gerstner waves and thei generalizations in hydrodynamics and geophysics, Phys.-Uspekhi, 65, 453-467.
- Bressan A, and Constantin A 2019 The Deflection Angle of Surface Ocean Currents From the Wind Direction, J. Geophys. Res.: Oceans, 124, 7412-7420.
- Chelton D B, deSzoeke R A, Schlax M G, El Naggar K, and Siwertz N 1998 Geographical Variability of the First Baroclinic Rossby Radius of Deformation, J. Phys. Oceanogr., 48, 433-460.
- Constantin A 2011 Nonlinear water waves with applications to wave-current interactions and tsunamis, Philadelphia: Society for Industrial and Applied Mathematics.
- Constantin A 2021 Frictional effects in wind-driven ocean currents, Geophys. Astrophys. Fluid Dyn., 115, 1–14.
- Constantin A 2022a Nonlinear wind-drift ocean currents in arctic regions, Geophys. Astrophys. Fluid Dyn., 116, 101–115.
- Constantin A 2022b Comments on: nonlinear wind-drift ocean currents in arctic regions, Geophys. Astrophys. Fluid Dyn., 116, 116–121.
- Constantin A, Dritschel D G, and Paldon N 2020 The deflection angle between a wind-forced surface current and the overlying wind in an ocean with vertically varying eddy viscosity, Phys. Fluids, 116604.
- Constantin A, and Ivanov R I 2019 Equatorial wave-current interactions, Comm. Math. Phys., 370, 1–48.
- Constantin A, and Johnson R S 2015 The dynamics of waves interacting with the equatorial undercurrent, Geophys. Astrophys. Fluid Dyn., 109, 311–358.
- Constantin A, and Johnson R S 2016 Current and future prospects for the application of systematic theoretical models to the study of problems in physical oceanography, Phys. Lett. A, 380, 3007–3012.
- Constantin A, and Johnson R S 2019a, Ekman-type solutions for shallow-water flows on a rotating sphere: A new perspective on a classical problem, Physics of Fluids, 31.2, 021401.
- Constantin A, and Johnson R S 2019b Large-scale oceanic currents as shallow-water asymptotic solutions of the Navier-Stokes equation in rotating spherical coordinates, Deep-Sea Res. II: Top. Stud. Oceanogr, 160, 32-40.
- Constantin A, and Johnson R S 2021, On the modelling of large-scale atmospheric flows, J. Differ. Equ., 285, 751–798.
- Constantin A, and Johnson R S 2023, On the dynamics of the near-surface currents in the Arctic Ocean, Nonlinear Anal.: Real World Appl. , 73, 103894.
- Constantin A, and Johnson R S 2024, Nonlinear Atmospheric Flow Patterns Confined to Zonal Cloud Bands, J. Atmos. Sci., 82.1, 3-15.
- Constantin A, and Monismith S G 2017, Gerstner waves in the presence of mean currents and rotation, J. Fluid Mech., 820, 511–528.
- Constantin A, and Strauss W A 2010, Trochoidal Solutions to the Incompressible Two-Dimensional Euler Equations, J. Math. Fluid Mech., 12, 181-201.
- Dritschel D G, Paldor N, and Constantin A 2020 The Ekman spiral for piecewise-uniform viscosity, Ocean Sci., 16, 1089-1093.
Ekman V W 1905 On the influence of the Earth's rotation on ocean-currents', Ark. Mat. Astr. Fys., 2, 1–52.
Henry D 2008 On Gerstner's Water Wave, J. Nonlinear Math. Phys., 15, 87-95.
Johnson R S 2017 Application of the ideas and techniques of classical fluid mechanics to some problems in physical oceanography, Phil. Trans. R. Soc. A, 376, 20170092.
Johnson R S 2019 Some problems in physical oceanography (including the use of rotating spherical coordinates) treated as exercises in classical fluid mechanics: Methods and examples, Deep-Sea Res. II: Top. Stud. Oceanogr, 160, 48-57.
Johnson R S 2022 The ocean and the atmosphere: An applied mathematician's view, Commun. Pure Appl. Math., 21 (7), 2357-2381.
Johnson R S 2023 An Introduction to the Mathematical Fluid Dynamics of Oceanic and Atmospheric Flows, Berlin: EMS Press.
Kamenkovich V M 1977 Fundamentals of Ocean Dynamics, Amsterdam: Elsevier.
Kundu P K 1990 Fluid Mechanics, San Diego: Academic Press.
Marshall J, and Plumb R 2016 A Atmosphere, Ocean and Climate Dynamics: An Introductory Text, New York: Academic Press.
Milne-Thomson L M 1938 Theoretical Hydrodynamics, London: The Macmillan Co..
Polton J A, Lenn Y-D, Elipot S, Chereskin T K, and Sprintall J 2013 Can Drake Passage observations match Ekman's classic theory?, J. Phys. Oceanogr., 43, 1733–1740.
Puntini C 2025 On the modeling of nonlinear wind-induced ice-drift ocean currents at the North Pole, arXiv preprint, arXiv:2503.12906.
Roach C J, Phillips H E, Bindoff N L, and Rintoul, S R 2015 Detecting and characterizing Ekman currents in the Southern Ocean, J. Phys. Oceanogr., 45, 1205–1223.
Roberti L 2021 Perturbation analysis for the surface deflection angle of Ekman-type flows with variable eddy viscosity, J. Math. Fluid Mech., 23, No. 57.
Roberti L 2022 The surface current of Ekman flows with time-dependent eddy viscosity, Comm. Pure Appl. Anal., 21, 2463–2477.
Roberti L 2025 Mathematical aspects of geophysical fluid flows, Ph.D. Thesis, University of Vienna, personal communication.
Smith S D 1988 Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature, J. Geophys. Res., 93, 15467–15472.
Soloviev A, and Lukas R 2014 The Near-Surface Layer of the Ocean: Structure, Dynamics and Applications, 2nd edn, Dordrecht: Springer.
Talley L D , Pickard G L, Emery W J, and Swift J H 2011 Descriptive Physical Oceanography: An Introduction, Amsterdam: Elsevier.
Trenberth J E, Large W G, Olson J G 1990 The mean annual cycle in global ocean wind stress, J. Phys. Oceanogr., 20, 1742-1760.
Vallis G K 2017 Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation, Cambridge: Cambridge University Press.
White A A 2002 A view of the equations of meteorological dynamics and various approximations, in J. Norbury & I. Roulstone (eds), Large-Scale Atmosphere-Ocean Dynamics, Cambridge: Cambridge University Press, 1–100.
Yoshikawa Y, and Masuda A 2009 Seasonal variations in the speed factor and deflection angle of the wind driven surface flow in the Tsushima Strait, J. Geophys. Res., 114, C12022.