β
1. Introduction
One of the classic problems in physical oceanography concerns the effects of wind forcing on the ocean. Its roots trace back to the pioneering work of Ekman (Reference Ekman1905), who identified the crucial role of the Earth’s rotation in explaining the observation made by Nansen during the Fram Arctic expedition of 1893–1896 (see Nansen Reference Nansen1897 for his own account of the voyage): steady wind-driven ocean currents are deflected to the right of the prevailing wind direction in the Northern Hemisphere. This deflection arises from the momentum balance between the Coriolis acceleration acting on the upper ocean layers and the frictional forces generated by the turbulent stress exerted by the wind.
The study of wind-driven flows has also played a fundamental role in the development of modern theories of large-scale ocean circulation. In particular, the seminal work of Sverdrup (Reference Sverdrup1947) established the relation between wind-stress curl and meridional mass transport in the interior of the ocean. Building on earlier ideas of Rossby (Reference Rossby1936) concerning lateral eddy viscosity, Munk (Reference Munk1950) derived a solution for the wind-driven circulation in a closed ocean basin and provided a theoretical explanation for the structure of large-scale ocean gyres. Together with the theory of western boundary intensification introduced by Stommel (Reference Stommel1948), these developments laid the foundations of dynamical oceanography and remain a cornerstone of our current understanding of the global ocean circulation. For a comprehensive account of these classical theories and their subsequent developments, we refer the reader to Pedlosky (Reference Pedlosky1996).
Despite the central importance of these ideas, theoretical investigations of wind-driven flows have become progressively less prominent in recent years. Increasing attention has been devoted to data-driven analyses and high-resolution numerical simulations, while analytical approaches based on classical fluid mechanics have received comparatively less emphasis (see, e.g. Chereskin Reference Chereskin1995; Niiler & Paduan Reference Niiler and Paduan1995; Poulain et al. Reference Poulain, Gerin, Mauri and Pennel2009; Yoshikawa & Masuda Reference Yoshikawa and Masuda2009; Yoshikawa et al. Reference Yoshikawa, Masuda, Marubayashi and Ishibashi2010; Kim et al. Reference Kim, Cornuelle and Terrill2010; Röhrs & Christensen Reference Röhrs and Christensen2015;Sentchev et al. Reference Sentchev, Forget and Fraunié2017; Berta et al. Reference Berta2018 for recent data-driven studies of Ekman flows). Needless to say, observational and numerical data are indispensable for testing hypotheses and validating theoretical predictions. However, we contend that a careful and systematic (asymptotic) analysis of the governing equations remains crucial for identifying the fundamental balances underlying the observed dynamics and for gaining a deeper understanding of the mechanisms driving such complex phenomena.
Furthermore, most theoretical studies of Ekman flows are formulated within the
$f$
-plane or
$\beta$
-plane approximations; see, for instance, Grisogono (Reference Grisogono1995), Lewis & Belcher (Reference Lewis and Belcher2004), Constantin et al. (Reference Constantin, Dritschel and Paldor2020), Constantin (Reference Constantin2021); Roberti (Reference Roberti2021, Reference Roberti2022), Marynets (Reference Marynets2022), Stefanescu (Reference Stefanescu2024) and Puntini (Reference Puntini2025, Reference Puntini2026). These approximations are obtained through the ad hoc neglect of selected terms in the governing equations and, although they have proved useful in many contexts, their range of validity is intrinsically limited, since the underlying assumptions are appropriate only for flows of sufficiently small horizontal extent. Therefore, it is desirable to have a model in full spherical coordinates that may be able to capture large-scale flows as well. Unfortunately, the full Navier–Stokes equations in spherical coordinates are intractable analytically; thus, it is necessary to make some approximations while trying to minimise the assumptions required. In this paper, we develop such a framework and apply it to the classical problem of wind-driven ocean currents.
The approach we will pursue relies on asymptotic expansions and has had a well-documented success in fluid mechanics; see, e.g. Balmforth & Liu (Reference Balmforth and Liu2004), Fernández-Nietoa et al. (Reference Fernández-Nietoa, Noble and Vila2010) and Chambon et al. (Reference Chambon, Freydier, Naaim and Vila2020) for applications to non-Newtonian fluids or Peregrine (Reference Peregrine1967), Johnson (Reference Johnson1980), Constantin & Johnson (Reference Constantin and Johnson2008) and Yuan et al. (Reference Yuan, Zhang, Wang, Chen and Johnson2025) in the context of water waves. In geophysical fluid dynamics, this method was used by Constantin & Johnson (Reference Constantin and Johnson2017) for the description of ocean gyres and has subsequently been successfully applied to a variety of problems; see Constantin & Johnson (Reference Constantin and Johnson2018, Reference Constantin and Johnson2019a , Reference Constantin and Johnsonb , Reference Constantin and Johnson2021, Reference Constantin and Johnson2023, Reference Constantin and Johnson2024a ) and Zhang & Zhang (Reference Zhang and Zhang2025). For a detailed review of similar asymptotic methods in oceanography, we refer the reader to Johnson (Reference Johnson2018, Reference Johnson2022).
In the present work, we proceed as follows. We begin with the Navier–Stokes equations, expressed in a rotating system of spherical coordinates, along with the equation of mass conservation and appropriate boundary conditions. Two assumptions are made a priori: we model turbulent stresses using an eddy-viscosity parametrisation, with the eddy viscosity assumed to depend only on depth, and we assume the density to vary only with depth too. The relevance of such depth-dependent viscosity profiles has been emphasised by Cronin & Kessler (Reference Cronin and Kessler2009) and Cronin & Tozuka (Reference Cronin and Tozuka2016). Then, we carefully rescale all the quantities appearing in the equations and identify two small parameters, which are independent of each other: the thin-shell parameter (
$\varepsilon$
) and the Rossby number (
$\mathscr{R}$
). These two parameters emerge naturally from the Navier–Stokes equations in spherical coordinates via the non-dimensionalisation. The former represents the ratio of a depth scale – in our case, the so-called Ekman depth – and the radius of the Earth, and is guaranteed to be small regardless of the nature of the flow; the latter is a well-known parameter in physical oceanography that is typically small for large-scale flows strongly affected by the Coriolis force (cf. the discussion by Kantha & Clayson Reference Kantha and Clayson2000). As noted by Constantin & Johnson (Reference Constantin and Johnson2019b
), performing an asymptotic expansion first in the thin-shell parameter and then in the Rossby number effectively separates the dynamical effects of the geometry (via the thin-shell parameter) and the properties of the flow due to the dominant time and length scales (via the Rossby number). A further advantage of the asymptotic approach adopted here is that it provides a systematic procedure for computing the governing equations at successive orders of approximation. In the present work, we derive the leading-order equations and obtain several explicit solutions, together with bounds for first-order correction with respect to
$\mathscr{R}$
. More importantly, the asymptotic hierarchy can, in principle, be continued to arbitrarily high order. Since the equations arising at each order are explicitly computable, numerical approximations of any desired accuracy may be obtained by retaining sufficiently many terms in the expansion.
We believe that the combination of mathematical rigour and physical consistency underlying our derivation and analysis provides a valuable contribution to the modelling of wind-driven ocean currents. In addition, the framework developed here may serve as a foundation for future investigations of large-scale ocean circulation. Indeed, the use of spherical coordinates offers the possibility of overcoming some of the limitations inherent in the classical
$f$
-plane and
$\beta$
-plane theories, and may therefore lead to a more accurate description of basin-scale and global ocean dynamics.
1.1. Plan of the paper
The paper is organised as follows. In § 2, after some remarks on the frame of reference and the coordinate system, we first present the governing equations and discuss in detail the boundary conditions, then proceed to rescaling all quantities and non-dimensionalising the equations and the boundary conditions, thereby identifying the thin-shell parameter and the Rossby number. Then, § 3 is devoted to the derivation of the leading-order governing equation through asymptotics: first, we make an asymptotic expansion of all variables with respect to the thin-shell parameter; then, within the leading-order equations thus derived, we perform a second asymptotic expansion – now with respect to the Rossby number. In § 4, we show how our formulation in spherical coordinates can be used to generalise the classic theory of the ocean circulation in the
$\beta$
-plane due to Sverdrup (Reference Sverdrup1947) to large-scale flows. In § 5, we derive the leading-order problem for the Ekman flow away from the equator by splitting the solution into a geostrophic and an Ekman part. In particular, in § 5.1 we prove the existence and uniqueness of the solution to the problem for the leading-order Ekman flow and show that this solution describes a classical Ekman spiral and analyse the surface deflection angle (that is, the angle between the direction of the wind and the induce wind-drift current on the surface); then, in § 5.2 we discuss some concrete examples of explicit eddy-viscosity profiles. Here we can explicitly solve the equation, draw numerical plots that help to visualise the flow and calculate the surface deflection angle using the theory developed in § 5.1. The proofs of the results of § 5.1 and the calculations of § 5.2 are deferred to Appendices A and B, respectively. Appendix C provides the derivation of the governing equations of the first-order correction to the leading-order equations with respect to the Rossby number
$\mathscr{R}$
, along with some ideas on how to provide a priori bounds on this first-order correction only in terms of the leading-order solution. Finally, Appendix D includes some three-dimensional plots of the Ekman spirals of § 5.2.
2. Governing equations
In this section, we will first present the equations of motion and the boundary conditions, then we will introduce an appropriate scaling and non-dimensionalise the problem. In the following, we use primes to denote physical (dimensional) quantities; we will remove them when we rescale the variables.
First, let us discuss the frame of reference in which the problem is formulated. A suitable coordinate system is the system of right-handed spherical coordinates
$(\varphi , \theta , r')$
, where
$\theta \in [-({\pi }/{2}), {\pi }/{2} ]$
is the angle of latitude,
$\varphi \in [0,2\pi )$
is the angle of longitude and
$r'\in [0,\infty )$
is the radial distance from the origin; see figure 1 for a sketch of these spherical coordinates. The coordinate change from these spherical coordinates to the standard Cartesian coordinates
$(x',y',z')$
is given by
\begin{equation} \left \{\begin{aligned} x' &= r' \cos \theta \cos \varphi ,\\[0.2em] y' &= r' \cos \theta \sin \varphi ,\\[0.2em] z' &= r' \sin \theta , \end{aligned}\right .\quad \text{with inverse}\quad \left \{ \begin{aligned} &\varphi = \tan ^{-1}\left (\frac {y'}{x'}\right )\!,\\[0.2em] &\theta = \sin ^{-1}\left (\frac {z'}{\sqrt {x'^2+y'^2+z'^2}}\right )\!,\\[0.2em] & r'=\sqrt {x'^2+y'^2+z'^2}. \end{aligned}\right . \end{equation}
With each point for which
$r\gt 0$
and
$\theta \notin \{\pm ({\pi }/{2})\}$
we can associate the usual local orthonormal basis for
$\mathbb{R}^3$
of unit vectors
$(\boldsymbol{e}_{\varphi }, \boldsymbol{e}_{\theta }, \boldsymbol{e}_r)$
:
$\boldsymbol{e}_{\varphi }$
points eastward,
$\boldsymbol{e}_{\theta }$
northward and
$\boldsymbol{e}_r$
radially. The velocity field is denoted by
$\boldsymbol{u}'=u'\boldsymbol{e}_{\varphi } + v'\boldsymbol{e}_{\theta } + w'\boldsymbol{e}_r$
. Note that we must take care to exclude the poles, where the tangent vectors
$\boldsymbol{e}_{\varphi }$
and
$\boldsymbol{e}_{\theta }$
are not defined, as a consequence of the hairy ball theorem (Renteln Reference Renteln2013); flows that extend to the poles may be treated, for instance, by rotating the coordinate system and moving the singularities to the equator, as done by Constantin & Johnson (Reference Constantin and Johnson2024b
) and Puntini (Reference Puntini2025), but we will not pursue this in the present paper. The vectors
$(\boldsymbol{e}_{\varphi }, \boldsymbol{e}_{\theta }, \boldsymbol{e}_r)$
are related to the Cartesian unit basis vectors
$(\boldsymbol{e}_1, \boldsymbol{e}_2, \boldsymbol{e}_3)$
through
\begin{equation} \left \{ \begin{aligned} &\boldsymbol{e_{\varphi }}=-\sin \varphi \,\boldsymbol{e}_1+\cos \varphi \,\boldsymbol{e}_2,\\[0.2em] &\boldsymbol{e_{\theta }}=-\cos \varphi \sin \theta \,\boldsymbol{e}_1-\sin \varphi \sin \theta \,\boldsymbol{e}_2+\cos \theta \,\boldsymbol{e}_3,\\[0.2em] &\boldsymbol{e}_{r}=\cos \varphi \cos \theta \,\boldsymbol{e}_1+\sin \varphi \cos \theta \,\boldsymbol{e}_2+\sin \theta \,\boldsymbol{e}_3. \end{aligned}\right . \end{equation}
Moreover, it is crucial to point out that, denoting by
$\boldsymbol{\varOmega }'$
the Earth’s constant angular velocity, with magnitude
$\varOmega ' = |\boldsymbol{\varOmega }'| \approx 7.29 \times 10^{-5} \,\mathrm{rad\,s^{-1}}$
, we choose the reference system so that it rotates about the
$z'$
axis, which is aligned with
$\boldsymbol{\varOmega }'$
, at angular speed
$\varOmega '$
; as a result, the Earth appears stationary. Although this choice is advantageous in many ways, it has the drawback of not being inertial, which – as discussed, for example, by Pedlosky (Reference Pedlosky1987), Phillips (Reference Phillips1987) and Vallis (Reference Vallis2017) – will give rise to additional terms (corresponding to apparent forces) in the momentum equations, namely, the Coriolis force
$2\boldsymbol{\varOmega }'\times \boldsymbol{u}'$
and the centrifugal force
$\boldsymbol{\varOmega }' \times (\boldsymbol{\varOmega }' \times \boldsymbol{r}')$
, where
$\boldsymbol{r}'$
is the position vector.
The classical spherical coordinate system.
Figure 1. Long description
A diagram of the classical spherical coordinate system. The diagram features a sphere with a coordinate system centered at point O. The sphere is intersected by a plane, creating a circle. Key points on the sphere are labeled N and S, representing the north and south poles, respectively. Point P is marked on the surface of the sphere. The diagram includes radial distance r’, and angles θ and φ, which define the position of point P. The unit vectors e_r, e_θ, and e_φ are shown, indicating the directions of the spherical coordinates. The axes x’, y’, and z’ are also depicted, representing the Cartesian coordinate system.
Since, for most geophysical flows, the Earth’s oblateness is not relevant (see the discussions by White Reference White2002 and Constantin & Johnson Reference Constantin and Johnson2021), we would like to approximate the Earth’s surface as a sphere of radius
$R' \approx 6371\, \mathrm{km}$
. However, some care is needed in order to minimise the dynamical error due to this approximation. Let us first define the geopotential
$\varPhi '$
as the sum of the gravitational potential and the term
$-({1}/{2})\varOmega '^2\ell ^2$
, where
$\ell$
is the distance of
$\boldsymbol{r}'$
from the axis of rotation. Over time, Earth has developed an equatorial bulge to counteract the centrifugal force, making the geopotential force (or effective gravity)
$\boldsymbol{g}' = -\boldsymbol{\nabla }\varPhi '$
act everywhere nearly perpendicularly to the Earth’s surface; in other words, the Earth’s surface is essentially a level set of the geopotential
$\varPhi '$
. Therefore, we choose to treat the level sets of the geopotential
$\varPhi$
(which, in reality, would be rotation ellipsoids of tiny eccentricity) as if they were perfect spheres. Thus, in this approximation, the horizontal component of effective gravity vanishes, and we replace a potentially significant dynamic error with a minor geometric error; the geopotential
$\varPhi '$
then depends solely on the vertical coordinate, and, for many applications, we can approximate
$\varPhi '(r')=g'r'$
. Given that the Earth’s oblateness is minimal, the spherical approximation provides a highly accurate representation of the spheroid.
2.1. The equations of motion
We can now present the equations of motion. In addition to the velocity field
$\boldsymbol{u}'=u'\boldsymbol{e}_{\varphi } + v'\boldsymbol{e}_{\theta } + w'\boldsymbol{e}_r$
introduced previously, we denote by
$\rho '$
the density and by
$p'$
the pressure. We assume that the density varies only with depth, while still treating the fluid as incompressible, since a pressure change of approximately
$500\,\mathrm{kPa}$
leads to only a
$0.025\,\%$
variation in the density of water (see Maslowe Reference Maslowe1986 and Vallis Reference Vallis2017). Therefore, the mass conservation equation reduces to the incompressibility condition, which in spherical coordinates reads
\begin{equation} \frac {1}{r' \cos \theta } \frac {\partial u'}{\partial \varphi } + \frac {1}{r' \cos \theta } \frac {\partial }{\partial \theta } \left ( v' \cos \theta \right ) + \frac {1}{r'^2} \frac {\partial }{\partial r'} \bigl ( r'^2 w' \bigr )=0. \end{equation}
In addition to incompressibility, in the equations of motions we assume the fluid flow to be laminar, with the turbulent effects parametrised by the dynamic (effective) horizontal and vertical viscosities
$\mu '_{{h}}$
and
$\mu '_{{v}}$
, which, like the density, are assumed to vary only with respect to the depth. We note that the concept of eddy viscosity represents a simple model for the closure of the turbulence problem, based on the hypothesis that Reynolds stresses are proportional to the mean velocity gradient (and thus depend on flow conditions) but do not represent a physical characteristic of the fluid (Pedlosky Reference Pedlosky1987; Pope Reference Pope2000); we will make additional remarks on the eddy viscosity in §§ 2.3 and 6. With these assumptions, the Navier–Stokes equations in spherical rotating coordinates are (see the derivation by Puntini Reference Puntini2025)
\begin{equation} \begin{aligned} &\left (\frac {\partial }{\partial t'}+\frac {u'}{r'\cos \theta }\frac {\partial }{\partial \varphi }+\frac {v'}{r' }\frac {\partial }{\partial \theta }+w'\frac {\partial }{\partial r'}\right )\begin{pmatrix} u' \\[0.2em] v' \\[0.2em] w' \end{pmatrix} \\[3pt] &\quad +\frac {1}{r'}\begin{pmatrix}-u'v'\tan \theta + u'w' \\[0.2em] u'^2\tan \theta +v'w'\\[0.2em] -u'^2-v'^2 \end{pmatrix}+2\varOmega '\begin{pmatrix} -v'\sin \theta +w'\cos \theta \\[0.2em] u'\sin \theta \\[0.2em] -u'\cos \theta \end{pmatrix}\\[3pt] &\quad = -\frac {1}{\rho '}\begin{pmatrix} \dfrac {1}{r'\cos \theta }\dfrac {\partial p'}{\partial \varphi }\\[15pt] \dfrac {1}{r'}\dfrac {\partial p'}{\partial \theta }\\[15pt] \dfrac {\partial p'}{\partial r} \end{pmatrix} -\frac {g'R^{\prime 2}}{r^{\prime 2}} \begin{pmatrix} 0\\[0.2em] 0\\[0.2em] 1 \end{pmatrix} + \frac {\mu '_{{v}}}{\rho '} \left (\frac {\partial ^2}{\partial r'^2}+ \frac {2}{r'}\frac {\partial }{\partial r'}\right )\begin{pmatrix} u'\\[0.2em] v'\\[0.2em] w' \end{pmatrix}\\[3pt] &\quad \quad + \frac {\mu '_{{h}}}{\rho ' r'^2}\left (\frac {1}{ \cos ^2 \theta } \frac {\partial ^2 }{\partial \varphi ^2}+\frac {\partial ^2}{\partial \theta ^2} - \tan \theta \frac {\partial }{\partial \theta }\right ) \begin{pmatrix} u'\\[0.2em] v'\\[0.2em] w' \end{pmatrix}\\[3pt] &\quad \quad -\frac {1}{\rho 'r'^2\cos ^2\theta }\begin{pmatrix} \mu '_{{h}} u\\[0.3em] \mu '_{{h}} v\\[0.3em] 2\mu '_{{v}}(w'\cos ^2\theta -v'\sin \theta \cos \theta ) \end{pmatrix}+\frac {2\mu '_{{h}}}{\rho 'r'^2}\frac {\partial }{\partial \theta }\begin{pmatrix} 0\\[0.2em] w'\\[0.2em] -v' \end{pmatrix}\\[3pt] &\quad \quad +\frac {2\mu '_{{h}}}{\rho 'r'^2\cos \theta }\frac {\partial }{\partial \varphi }\begin{pmatrix} w'-v'\tan \theta \\[0.2em] u'\tan \theta \\[0.2em] -u' \end{pmatrix}+\frac {1}{\rho '} \frac {\mathrm{d} \mu '_{{v}}}{\mathrm{d} r'} \begin{pmatrix} r' \dfrac {\partial }{\partial r'} \left(\dfrac {u'}{r'}\right)\\[15pt] r' \dfrac {\partial }{\partial r'} \left(\dfrac {v'}{r'} \right)\\[15pt] 2 \dfrac {\partial w'}{\partial r'} \end{pmatrix}+\frac {1}{\rho '} \frac {\mathrm{d}\mu '_{{h}}}{\mathrm{d} r'} \begin{pmatrix} \dfrac {1}{r'\cos \theta }\dfrac {\partial {w'}}{\partial \varphi }\\[15pt] \dfrac {1}{r'} \dfrac {\partial {w'}}{\partial \theta }\\[15pt] {0} \end{pmatrix}\!, \end{aligned} \end{equation}
where
$R' \approx 6371\, \mathrm{km}$
is the Earth’s average radius and
$g' \approx 9.81 \, \mathrm{m\, s^{-2}}$
is the average gravitational acceleration at ground level. It is worth remarking that the Navier–Stokes equation (2.4), derived by Puntini (Reference Puntini2025) to correct some errors generally present in the literature, could be further generalised. Indeed, the tensor related to the eddy viscosity considered by Puntini (Reference Puntini2025) is isotropic in the horizontal directions and is assumed to vary only with depth. Accounting for a more general eddy-viscosity tensor (in the vein of recent works in the
$f$
-plane approximation, such as those by Wirth (Reference Wirth2010) or Bykov & Gordin Reference Bykov and Gordin2024), while technically possible and undoubtedly interesting, requires extensive calculations that would exceed the scope of this work. Therefore, for the time being, we content ourselves with the simpler stress tensor from Puntini (Reference Puntini2025) and leave the more general case for future investigations.
2.2. Boundary conditions
Let us now discuss the boundary conditions that characterise the dynamics of wind-driven ocean currents. At the free surface of the ocean, which we assume to be described by the equation
\begin{equation} r'=R'+h'(\varphi ,\theta , t') \end{equation}
for an unknown function
$h'$
, we impose the dynamic boundary condition
\begin{equation} p'=p_{{s}}\quad \text{on} \quad \{r'=R'+h'\}, \end{equation}
where
$p_{{s}}$
is the given (and, in general, non-constant) surface (atmospheric) pressure, and the kinematic boundary condition
\begin{equation} w'=\frac {\partial h'}{\partial t'}+\frac {u'}{r'\cos \theta }\frac {\partial h'}{\partial \varphi }+\frac {v'}{r'}\frac {\partial h'}{\partial \theta }\quad \text{on} \quad \{r'=R'+h'\}, \end{equation}
which ensures that no fluid particles escape the fluid via the free surface. On the rigid bottom, which for simplicity we assume to be flat, we require the velocity field to vanish
\begin{equation} \boldsymbol{u}' = 0 \quad \text{on} \quad \{r'=R'+H'\}, \end{equation}
for a constant
$H' \lt h'$
.
The stress generated by the wind blowing over the ocean’s surface
$\{r' = R'+h'\}$
is expressed by the bulk formula (Soloviev & Lucas Reference Soloviev and Lucas2014)
\begin{equation} \boldsymbol{\tau }'_{{a}}=\rho '_{{a}} C_{{aw}} |\boldsymbol{u}'_{{a}}-\boldsymbol{u}'|(\boldsymbol{u}'_{{a}}-\boldsymbol{u}') \quad \text{on} \quad \{r'=R'+h'\}, \end{equation}
where
$\rho '_{{a}} \approx 1.2\,$
kg m
$^{-3}$
is the density of the air,
$C_{{aw}}\approx 10^{-3}$
is the non-dimensional drag coefficient of air over water and
$\boldsymbol{u}'_{{a}}= u'_{{a}}\boldsymbol{e}_{\varphi } + v'_{{a}}\boldsymbol{e}_{\theta }$
is the horizontal wind velocity, conventionally measured at
$10\,\mathrm{m}$
above the sea surface. See figure 2 and Trenberth et al. (Reference Trenberth, Large and Olson1990) for data on the mean wind velocity over the ocean.
Monthly averaged ocean wind speed and direction vectors, with vector lengths proportional to the reference scale (in
$\mathrm{m\, s^{-1}}$
), based on observations from NASA’s QuikSCAT satellite. Image credit: NOAA.
The bulk formula (2.9) provides us with an additional boundary condition for our problem. In fact, on the surface of the ocean, the shear stress is given by
\begin{equation} \boldsymbol{\tau }'_{{w}}=\mu '_{{v}}\left (\frac {\partial u'}{\partial r'}\boldsymbol{e}_{\varphi } + \frac {\partial v'}{\partial r'}\boldsymbol{e}_{\theta }\right ) \quad \text{on} \quad \{r'=R'+h'\}, \end{equation}
where
$\mu '_{{v}}$
is the vertical dynamic viscosity, and consequently, in view of (2.9), the continuity of the stress across the water surface yields
\begin{equation} \mu '_{{v}}\begin{pmatrix} \dfrac {\partial u'}{\partial r'}\\[1em] \dfrac {\partial v'}{\partial r'} \end{pmatrix} =\rho '_{{a}} C_{{aw}} \sqrt {(u_{{a}}'-u')^2+(v_{{a}}'-v')^2}\begin{pmatrix} u'_{{a}}-u'\\[0.2em] v'_{{a}}-v' \end{pmatrix} \quad \text{on} \quad \{r'=R'+h'\}. \end{equation}
2.3. Scaling and non-dimensionalisation
Clearly, the Navier–Stokes equations (2.4) are intractable by analytical means in their full generality. In order to simplify them, we will be guided by the scaling of the variables appearing in the equations; through this, we will be able to identify the terms in the equations that contribute the most and those that, on the contrary, may be neglected with little error. To this end, we must identify the relevant scales of the flow.
Before we proceed, some clarifications on the concept of eddy viscosity are necessary, because the physical interpretation of a local eddy viscosity at the ocean surface is not straightforward. Near the air–sea interface, wave breaking and other surface processes can induce non-local momentum transport, so that turbulent stresses are not necessarily related to the local mean shear through a simple diffusive closure (Craig & Banner Reference Craig and Banner1994). Furthermore, the specification of an eddy viscosity at the ocean surface is inherently ambiguous. As an example, K-profile parametrisation (KPP) turbulence closures predict that the mixing length, and thus the eddy viscosity, vanishes as one approaches the surface (Large et al. Reference Large, McWilliams and Doney1994), whereas field observations suggest that wave-induced processes may substantially enhance near-surface mixing beyond the predictions of such models (Fisher et al. Reference Fisher, Sanford and Scully2018). At the same time, direct observations in the immediate vicinity of the air–sea interface are extremely challenging. Consequently, no single universally accepted value can be assigned to the surface eddy viscosity, and estimates may vary depending on the modelling framework and the representation of near-surface processes. In this work, the adopted value should therefore be interpreted as a convenient reference scale for non-dimensionalisation within the present modelling framework, rather than as a directly measurable physical quantity.
With this caveat, let us denote by
$\bar {\mu }'_{{h}}$
and
$\bar {\mu }'_{{v}}$
the reference values for the horizontal and vertical dynamic eddy viscosity on the surface, respectively. A typical value for
$\bar {\mu }'_{{v}}$
in the near-surface ocean is
$\bar {\mu }'_{{v}} \approx 10\,\mathrm{kg\,m^{-1}\,s^{-1}}$
(Kamenkovich Reference Kamenkovich1977; Wenegrat et al. Reference Wenegrat, McPhaden and Lien2014), although observations indicate substantial variability (Huang Reference Huang1979). A reference value for the horizontal (dynamic) eddy viscosity is even more difficult to determine; nevertheless, horizontal mixing greatly exceeds vertical mixing, with
$\bar {\mu }'_{{h}}/\bar {\mu }'_{{v}}\sim 10^3$
(Pedlosky Reference Pedlosky1987; Talley et al. Reference Talley, Pickard, Emery and Swift2011; Constantin & Johnson Reference Constantin and Johnson2018; Constantin Reference Constantin2022), implying
$\bar {\mu }'_{{h}}\gg \bar {\mu }'_{{v}}$
. However, because we non-dimensionalise the equations, their precise values are not essential. Furthermore, we denote by
\begin{equation} D'=\sqrt {\frac {\bar {\mu '_{{v}}} }{\bar {\rho }'\varOmega '}}\approx 12\, \mathrm{m} \end{equation}
the Ekman depth and by
$\bar {\rho }' = 1.025\,\mathrm{kg\,m^{-3}}$
the average seawater density at the surface. The choice for the speed scale
$U'$
arises naturally from (2.11) for the wind stress on the surface. In fact, using (2.16) and setting
$(u_{{a}}',v_{{a}}') = U_{{a}}'(u_{{a}},v_{{a}})$
for an average wind speed
$U_{{a}}$
(for instance,
$U_{{a}}' = 10\,\mathrm{m\,s^{-1}}$
), it follows immediately that the appropriate scale for
$U'$
is
\begin{equation} U' = \frac {\rho '_{{a}} C_{{aw}} U_{{a}}^{\prime 2}D'}{\bar {\mu }'_{{v}}}. \end{equation}
Substituting the values mentioned above and taking
$U'_{{a}} = 10\,\mathrm{m\,s^{-1}}$
, we find
\begin{equation} U' \approx 0.14\,\mathrm{m\,s^{-1}}. \end{equation}
Based on these scales, we introduce two non-dimensional parameters
\begin{equation} \varepsilon =\frac {D'}{R'}\approx 2 \times 10^{-6} \quad \text{and} \quad \mathscr{R} = \frac {U'}{\varOmega ' R'}\approx 2 \times 10^{-4}; \end{equation}
the former is the thin-shell parameter (the ratio of the Ekman length and the Earth’s average radius), and the latter is the Rossby number (a well-known parameter in oceanography, defined as the ratio between the scales of the advective terms and the Coriolis terms in the Navier–Stokes equations). Intuitively, a small Rossby number indicates that the advective terms are dominated by the Coriolis terms; this is usually the case for large-scale oceanic flows (cf. the discussion by Kantha & Clayson Reference Kantha and Clayson2000 and Vallis Reference Vallis2017).
In view of these considerations, we can non-dimensionalise the equations (2.3) and (2.4) by setting
\begin{equation} \begin{gathered} r'=R'+D'z, \quad (u',v',w')=U'(u,v,\varepsilon w), \quad \rho '(r') = \bar {\rho }'\rho (z), \\ p'= \bar {\rho }'U'\varOmega 'R'p, \quad t'=\frac {R'}{U'}\,t,\quad \mu '_{{h}}(r')=\bar {\mu }'_{{h}} n(z), \quad \mu '_{{v}}(r')=\bar {\mu }'_{{v}}m(z). \end{gathered} \end{equation}
The factor
$\varepsilon$
in front of
$w$
ensures the consistency of the scaling. Indeed, given the particle trajectories
$t'\mapsto \varPhi (t') = (\varphi (t'),\theta (t'),r'(t'))$
in spherical coordinates, we have
$w'(\varPhi (t'),t') = {\mathrm{d} r'(t')}/{\mathrm{d} t'}$
, hence, taking (2.16) into account and setting
$w = {\mathrm{d} z}/{\mathrm{d} t}$
, it follows
$w' = \varepsilon U'w$
.
With (2.16), the first two components of (2.4) are
\begin{align} &\mathscr{R}\left (\frac {\partial }{\partial t}+\frac {u}{(1+\varepsilon z)\cos \theta }\frac {\partial }{\partial \varphi }+\frac {v}{1+\varepsilon z}\frac {\partial }{\partial \theta }+w\frac {\partial }{\partial z}\right )\begin{pmatrix} u \\[0.2em] v \end{pmatrix} \nonumber\\[5pt] &+\frac {\mathscr{R}}{1+\varepsilon z}\begin{pmatrix}-uv\tan \theta + \varepsilon uw \\[0.2em] u^2\tan \theta +\varepsilon vw \end{pmatrix}+2 \begin{pmatrix} -v\sin \theta +\varepsilon w\cos \theta \\[0.2em] u\sin \theta \end{pmatrix} \nonumber\\[5pt] &\quad = -\frac {1}{1+\varepsilon z}\begin{pmatrix}\dfrac {1}{\rho \cos \theta }\dfrac {\partial p}{\partial \varphi } \\[15pt] \dfrac {1}{\rho }\dfrac {\partial p}{\partial \theta }\end{pmatrix}+ \frac {m}{\rho } \left (\frac {\partial ^2}{\partial z^2}+ \frac {2\varepsilon }{(1+\varepsilon z)}\frac {\partial }{\partial z}\right )\begin{pmatrix} u \\[0.2em] v \end{pmatrix} \nonumber\\[5pt] &\quad \quad +\frac {{\varepsilon ^2\unicode {x03BD}} n}{\rho (1+\varepsilon z)^2} \left (\frac {1}{ \cos ^2 \theta } \frac {\partial ^2 }{\partial \varphi ^2}+\frac {\partial ^2}{\partial \theta ^2} - \tan \theta \frac {\partial }{\partial \theta }\right ) \begin{pmatrix} u\\[0.2em] v \end{pmatrix} \nonumber\\[5pt] &\quad \quad -\frac {{\varepsilon ^2\unicode {x03BD}} n}{\rho (1+\varepsilon z)^2\cos ^2\theta }\begin{pmatrix} u\\[0.2em] v \end{pmatrix}+\frac {2 {\varepsilon ^2\unicode {x03BD}} n}{\rho (1+\varepsilon z)^2} \left (\frac {\partial }{\partial \theta }\begin{pmatrix} 0\\[0.2em] \varepsilon w \end{pmatrix}+\frac {1}{\cos \theta }\frac {\partial }{\partial \varphi }\begin{pmatrix} \varepsilon w-v\tan \theta \\[0.2em] u\tan \theta \end{pmatrix}\right ) \nonumber\\[5pt] &\quad \quad + \frac {1}{\rho }\frac {\mathrm{d} m}{\mathrm{d} z} \begin{pmatrix} \dfrac {\partial u}{\partial z}-\dfrac {u\varepsilon }{(1+\varepsilon z)^2}\\[15pt] \dfrac {\partial v}{\partial z}-\dfrac {v\varepsilon }{(1+\varepsilon z)^2} \end{pmatrix}+ \frac {\varepsilon ^{2}\unicode {x03BD}}{\rho (1+\varepsilon z)}\frac {\mathrm{d} n}{\mathrm{d} z} \begin{pmatrix} \dfrac {1}{\cos \theta }\dfrac {\partial {w}}{\partial \varphi }\\[15pt] \dfrac {\partial {w}}{\partial \theta } \end{pmatrix}\!, \end{align}
where we denoted the ratio between the scales of the dynamic eddy viscosities by
$\unicode {x03BD}= {\bar {\mu '_{{h}}} }/{\bar {\mu '_{{v}}} }$
, while, writing
\begin{equation} g = \frac {g'D'}{U'\varOmega 'R'} =\mathcal{O}(1), \end{equation}
the third component of (2.4) is
\begin{equation} \begin{aligned} &\varepsilon ^2\mathscr{R}\left (\frac {\partial w}{\partial t}+\frac {u}{(1+\varepsilon z)\cos \theta }\frac {\partial w}{\partial \varphi }+\frac {v}{(1+\varepsilon z)}\frac {\partial w}{\partial \theta }+ w\frac {\partial w}{\partial z}\right ) \nonumber \\& \quad \quad -\frac {\varepsilon \mathscr{R}}{1+\varepsilon z} (u^2+v^2)-2 u \varepsilon \cos \theta \\&\quad = - \frac {1}{\rho }\frac {\partial p}{\partial z} - \frac {g}{(1+\varepsilon z)^2} + \varepsilon ^2 \frac {m}{\rho } \left (\frac {\partial ^2 w}{\partial z^2}+ \frac {2\varepsilon }{(1+\varepsilon z)}\frac {\partial w}{\partial z}\right ) \\ &\quad \quad +\frac {{\varepsilon ^{4}\unicode {x03BD}}n}{\rho (1+\varepsilon z)^2}\left (\frac {1}{ \cos ^2 \theta } \frac {\partial ^2 w }{\partial \varphi ^2}+\frac {\partial ^2 w}{\partial \theta ^2} - \tan \theta \frac {\partial w}{\partial \theta }\right )\\&\quad \quad -\frac {2\varepsilon ^3m}{\rho (1+\varepsilon z)^2\cos ^2\theta }(\varepsilon w \cos ^2\theta -v\sin \theta \cos \theta ) \\ &\quad \quad -\frac {2 {\varepsilon ^3\unicode {x03BD}} n}{\rho (1+\varepsilon z)^2} \left (\frac {\partial v}{\partial \theta } +\frac {1}{\cos \theta }\frac {\partial u}{\partial \varphi } \right ) + 2\frac {\varepsilon ^2}{\rho } \frac {\mathrm{d} m}{\mathrm{d} z} \frac {\partial w}{\partial z}. \end{aligned} \end{equation}
Similarly, from the continuity equation (2.3) we end up with the non-dimensional equation
\begin{equation} (1+\varepsilon z)\left [ \frac {\partial u}{\partial \varphi } + \frac {\partial }{\partial \theta } \left ( v \cos \theta \right )\right ] + \cos \theta \frac {\partial }{\partial z} \bigl ( (1+\varepsilon z)^2 w \bigr )=0. \end{equation}
Analogously to the equations of motion, we have to non-dimensionalise the boundary conditions as well. In order to do so, we start from the kinematic boundary condition. Setting
$h'=D'h$
and
$H' = D'H$
, (2.7) reads
\begin{equation} w = \frac {\partial h}{\partial t}+\frac {u}{(1+\varepsilon z)\cos \theta }\frac {\partial h}{\partial \varphi }+\frac {v}{(1+\varepsilon z)}\frac {\partial h}{\partial \theta }\quad \text{on} \quad \{z=h\}, \end{equation}
while, setting
$\boldsymbol{u} = u\boldsymbol{e}_{\varphi } + v\boldsymbol{e}_{\theta } + w\boldsymbol{e}_r$
, (2.8) becomes
\begin{equation} \boldsymbol{u} = 0 \quad \text{on} \quad \{z=H\}. \end{equation}
For the dynamic boundary condition, in view of (2.6) and (2.16), on the free surface we have
\begin{equation} p = p_{{s}} \quad \text{on} \quad \{z=h\}, \end{equation}
where
\begin{equation} p_{{s}} = \frac {p'_{{s}}}{\rho 'U'\varOmega 'R'}. \end{equation}
Finally, setting
\begin{equation} \sigma = \frac {U'}{U_{{a}}'} \approx 0.014, \end{equation}
from (2.13) we see that the non-dimensional form of (2.11) is
\begin{equation} m\begin{pmatrix} \dfrac {\partial u}{\partial z}\\[1em] \dfrac {\partial v}{\partial z} \end{pmatrix} = \sqrt {(u_{{a}}-\sigma u)^2+(v_{{a}}-\sigma v')^2}\begin{pmatrix} u_{{a}}-\sigma u\\[0.2em] v_{{a}}-\sigma v \end{pmatrix} \quad \text{on} \quad \{z=h\}. \end{equation}
3. Derivation of the leading-order equations
Now that the equations have been non-dimensionalised, we would like to exploit the smallness of the parameters
$\varepsilon$
and
$\mathscr{R}$
to simplify the problem using asymptotic expansions, which should provide us at least with viable approximations of the full solution; for more details on the general theory, we refer the reader to Van Dyke (Reference Van Dyke1975), Johnson (Reference Johnson2005) and Holmes (Reference Holmes2015). First, let us forget about
$\mathscr{R}$
for the moment and make the ansatz
\begin{equation} q(\varphi ,\theta ,z,t;\varepsilon ) = q_0(\varphi ,\theta ,z,t) + o(\varepsilon ) \quad \text{as} \quad \varepsilon \searrow 0, \end{equation}
where
$q$
(and, correspondingly,
$q_0$
) represents each of the variables
$u,\,v,\,w,\,P,\,h$
. This ansatz is motivated by the structure of (2.17) and (2.19), where
$\varepsilon$
appears only in integer powers, as well as through
$(1+\varepsilon z)^{-1} = 1 - \varepsilon z + o(\varepsilon )$
. Plugging (3.1) into (2.17) and equating the coefficients of the various powers of
$\varepsilon$
, we obtain, at leading order with respect to
$\varepsilon$
,
\begin{equation} \begin{aligned} &\mathscr{R}\left \{\left (\frac {\partial }{\partial t}+\frac {u_0}{\cos \theta }\frac {\partial }{\partial \varphi }+v_0\frac {\partial }{\partial \theta }+w_0\frac {\partial }{\partial z}\right )\begin{pmatrix} u_0 \\[0.2em] v_0 \end{pmatrix} + \begin{pmatrix}-u_0v_0\tan \theta \\[0.2em] u_0^2\tan \theta \end{pmatrix}\right \} \\ &+ 2 \begin{pmatrix} -v_0\sin \theta \\[0.2em] u_0\sin \theta \end{pmatrix} + \begin{pmatrix}\dfrac {1}{\rho \cos \theta }\dfrac {\partial p_0}{\partial \varphi } \\[15pt] \dfrac {1}{\rho }\dfrac {\partial p_0}{\partial \theta }\end{pmatrix} - \frac {1}{\rho }\frac {\partial }{\partial z}\left (m\frac {\partial }{\partial z}\begin{pmatrix} u_0 \\[0.2em] v_0 \end{pmatrix}\right ) = 0, \end{aligned} \end{equation}
while from (2.19) we get
\begin{equation} \frac {\partial p_0}{\partial z} = -g\rho, \end{equation}
and the continuity equation (2.20) yields
\begin{equation} \frac {\partial u_0}{\partial \varphi } + \frac {\partial }{\partial \theta } \left (v_0 \cos \theta \right ) + \cos \theta \frac {\partial w_0}{\partial z} = 0. \end{equation}
From the kinematic boundary conditions (2.21) we obtain
\begin{equation} w_0 = \frac {\partial h_0}{\partial t}+\frac {u_0}{\cos \theta }\frac {\partial h_0}{\partial \varphi }+v_0\frac {\partial h_0}{\partial \theta } \quad \text{on} \quad \{z=h_0\}. \end{equation}
For simplicity, henceforth we assume that
$h_0$
is a simple flat surface, that is,
\begin{equation} h_0\equiv 0, \end{equation}
which implies that
$w_0=0$
there. Writing
$\boldsymbol{u}_0 = u_0\boldsymbol{e}_{\varphi } + v_0\boldsymbol{e}_{\theta } + w_0\boldsymbol{e}_{\theta }$
, (2.22) becomes
\begin{equation} \boldsymbol{u}_0 = 0 \quad \text{on} \quad \{z=H\}. \end{equation}
With this assumption on the free surface, the dynamic boundary condition (2.23) is simply
\begin{equation} p_0 = p_{{s}} \quad \text{on } \quad \{z=0\}, \end{equation}
and the upper boundary condition (2.26) for the wind stress gives
\begin{equation} \begin{pmatrix} \dfrac {\partial u_0}{\partial z}\\[1em] \dfrac {\partial v_0}{\partial z} \end{pmatrix} = \sqrt {(u_{{a}}-\sigma u_0)^2+(v_{{a}}-\sigma v_0)^2}\begin{pmatrix} u_{{a}}-\sigma u_0\\[0.2em] v_{{a}}-\sigma v_0 \end{pmatrix}\quad \text{on} \quad \{z=0\}, \end{equation}
noting that, with our scaling,
$m=1$
on
$\{z=0\}$
.
The next step is to perform a second asymptotic expansion, this time with respect to
$\mathscr{R}$
, within the leading-order equations with respect to
$\varepsilon$
that we have derived. To this end, let us make the ansatz
\begin{equation} q_0\bigl (\varphi ,\theta ,z,t;\mathscr{R}\bigr )\sim \tilde {q}(\varphi ,\theta ,y,t) + \mathscr{R}\hat {q}(\varphi ,\theta ,y,t) + o(\mathscr{R}\,), \end{equation}
where, as in (3.1),
$q$
denotes each of the variables
$u,\,v,\,w,$
and
$P$
. Plugging (3.10) into (3.2) and collecting terms of the same order in
$\mathscr{R}$
, we see that
\begin{equation} \begin{aligned} &\left \{2 \begin{pmatrix} -\tilde {v}\sin \theta \\[0.2em] \tilde {u}\sin \theta \end{pmatrix} + \begin{pmatrix} \dfrac {1}{\rho \cos \theta }\dfrac {\partial \tilde {p}}{\partial \varphi }\\[15pt] \dfrac {1}{\rho }\dfrac {\partial \tilde {p}}{\partial \theta } \end{pmatrix} - \frac {1}{\rho }\frac {\partial }{\partial z}\left (m\frac {\partial }{\partial z}\begin{pmatrix} \tilde {u}\\[0.2em] \tilde {v} \end{pmatrix}\right )\right \} \\[3pt] &+ \mathscr{R}\left \{\left (\frac {\partial }{\partial t} + \frac {\tilde {u}}{\cos \theta }\frac {\partial }{\partial \varphi } + \tilde {v}\frac {\partial }{\partial \theta } + \tilde {w}\frac {\partial }{\partial z}\right )\begin{pmatrix} \tilde {u}\\[0.2em] \tilde {v} \end{pmatrix} + \begin{pmatrix} -\tilde {u}\tilde {v}\tan \theta \\[0.2em] \tilde {u}^2\tan \theta \end{pmatrix} \right .\\[3pt] & +\biggl . 2 \begin{pmatrix} -\hat {v}\sin \theta \\[0.2em] \hat {u}\sin \theta \end{pmatrix} + \begin{pmatrix} \dfrac {1}{\rho \cos \theta }\dfrac {\partial \hat {p}}{\partial \varphi }\\[15pt] \dfrac {1}{\rho }\dfrac {\partial \hat {p}}{\partial \theta } \end{pmatrix} - \frac {1}{\rho }\frac {\partial }{\partial z}\left (m\frac {\partial }{\partial z}\begin{pmatrix} \hat {u}\\[0.2em] \hat {v} \end{pmatrix}\right )\biggr \} + o(\mathscr{R}\,) = 0, \end{aligned} \end{equation}
while doing the same in (3.3) yields
\begin{equation} \frac {\partial \tilde {p}}{\partial z} + g\rho + \mathscr{R}\frac {\partial \hat {p}}{\partial z} + o(\mathscr{R}\,) = 0. \end{equation}
Performing the same procedure in the continuity equation (3.4) leads us to
\begin{equation} \frac {\partial \tilde {u}}{\partial \varphi } + \frac {\partial }{\partial \theta } \left ( \tilde {v} \cos \theta \right ) + \cos \theta \frac {\partial \tilde {w}}{\partial z} + \mathscr{R}\left \{\frac {\partial \hat {u}}{\partial \varphi } + \frac {\partial }{\partial \theta } \left (\hat {v} \cos \theta \right ) + \cos \theta \frac {\partial \tilde {w}}{\partial z}\right \} + o(\mathscr{R}\,) = 0, \end{equation}
whereas, for the dynamic boundary condition, we have
\begin{equation} \tilde {p} + \mathscr{R}\hat {p} + o(\mathscr{R}\,) = 0 \quad \text{on} \quad \{z=0\}. \end{equation}
For the wind stress, inserting (3.10) into (3.9), and using the Taylor expansion
$\sqrt {1+x} = 1 +( {1}/{2})x + o(x)$
as
$x\to 0$
, we see that
\begin{equation} \begin{aligned} &\sqrt {(u_{{a}}-\sigma u_0)^2+(v_{{a}}-\sigma v_0)^2} \\[3pt] &\hspace {2cm}= \sqrt {(u_{{a}}-\sigma \tilde {u})^2+(v_{{a}}-\sigma \tilde {v})^2} \\[3pt] &\hspace {2cm} - \mathscr{R}\sigma \frac {\hat {u}(u_{{a}} - \sigma \tilde {u}) + \hat {v}(v_{{a}} - \sigma \tilde {v})}{\sqrt {(u_{{a}}-\sigma \tilde {u})^2+(v_{{ a}}-\sigma \tilde {v})^2}} + o(\mathscr{R}\,) \quad \text{on} \quad \{z=0\}. \end{aligned} \end{equation}
Setting coefficients to zero in (3.11), (3.12), and (3.13), we obtain the leading-order set of equations
\begin{equation} \left .\begin{aligned} & \frac {1}{\rho }\frac {\partial }{\partial z}\left (m\frac {\partial \tilde {u}}{\partial z}\right ) + 2\tilde {v}\sin \theta = \frac {1}{\rho \cos \theta }\frac {\partial \tilde {p}}{\partial \varphi } \\[11pt] &\frac {1}{\rho }\frac {\partial }{\partial z}\left (m\frac {\partial \tilde {v}}{\partial z}\right ) - 2\tilde {u}\sin \theta = \frac {1}{\rho }\frac {\partial \tilde {p}}{\partial \theta } \\[11pt] & \frac {\partial \tilde {p}}{\partial z} = -g\rho \\[11pt] &\frac {\partial \tilde {u}}{\partial \varphi } + \frac {\partial }{\partial \theta } \left ( \tilde {v} \cos \theta \right ) + \cos \theta \frac {\partial \tilde {w}}{\partial z} = 0 \end{aligned}\;\right \} \quad \text{in} \quad \{H \lt z \lt 0\}; \end{equation}
recall that here
$m$
denotes the non-dimensional eddy viscosity. The condition (3.7) reads simply
\begin{equation} \tilde {u}\boldsymbol{e}_{\varphi } + \tilde {v}\boldsymbol{e}_{\theta } + \tilde {w}\boldsymbol{e}_r = 0 \quad \text{on} \quad \{z=H\}; \end{equation}
the dynamic boundary condition (3.14) with (2.23) yields
\begin{equation} \tilde {p} = p_{{s}} \quad \text{on} \quad \{z=0\}; \end{equation}
and, from (3.15), we obtain at leading order
\begin{equation} \begin{pmatrix} \dfrac {\partial \tilde {u}}{\partial z}\\[12pt] \dfrac {\partial \tilde {v}}{\partial z} \end{pmatrix} = \sqrt {(u_{{a}}-\sigma \tilde {u})^2+(v_{{a}}-\sigma \tilde {v})^2}\begin{pmatrix} u_{{a}}-\sigma \tilde {u}\\[0.2em] v_{{a}}-\sigma \tilde {v} \end{pmatrix} =: \begin{pmatrix} \tau _\varphi \\[0.2em] \tau _\theta \end{pmatrix} \quad \text{on} \quad \{z=0\}. \end{equation}
Finally, (3.7) yields
\begin{equation} \tilde {u}\boldsymbol{e}_{\varphi } + \tilde {v}\boldsymbol{e}_{\theta } + \tilde {w}\boldsymbol{e}_r = 0 \quad \text{on} \quad \{z=H\}. \end{equation}
The first two equations in (3.16) are essentially a second-order system of ordinary differential equations (ODEs) with respect to the variable
$z$
, in which the other variables
$\varphi$
and
$\theta$
appear only parametrically. Physically, they describe a three-way balance between the Coriolis force, the horizontal pressure gradient force and the viscous stress caused by the wind.
4. Sverdrup’s theory of the ocean circulation: Comparison with
$\beta$
-plane theory
Before presenting results on the Ekman currents, we briefly digress to consider another problem for which the present asymptotic expansion of the equations of motion in spherical coordinates may prove particularly useful: the understanding of large-scale ocean circulation. An extensive application of the present framework to this problem is ongoing and will be the subject of future works. Nevertheless, to provide some insight into its potential, we consider the simplest model of ocean circulation, namely that of Sverdrup (Reference Sverdrup1947). Even at the level of this idealised model, the present equations – which retain the spherical geometry of the Earth at leading order by allowing both
$\theta$
and
$\varphi$
to vary – lead to results that differ in certain respects from the classical picture proposed by Sverdrup (Reference Sverdrup1947).
Following Sverdrup (Reference Sverdrup1947) (see also Stewart Reference Stewart2009), we define
$\varPi$
,
$M_1$
and
$M_2$
by
\begin{align} \frac {\partial \varPi }{\partial \varphi } = \int ^0_H\frac {\partial \tilde {p}}{\partial \varphi }\,\textrm {d}z, \quad \frac {\partial \varPi }{\partial \theta } = \int ^0_H\frac {\partial \tilde {p}}{\partial \theta }\,\textrm {d}z,\quad M_1=\int _H^0 \rho \tilde {u}\,\textrm {d}z \quad \text{and} \quad M_2=\int _H^0 \rho \tilde {v}\,\textrm {d}z, \end{align}
respectively, where, as before,
$H$
is the bottom of the ocean (corresponding, in dimensional variables, to some thousands of metres). Following Sverdrup, we may argue that at such depth
$H$
the stress is negligible, as the total current goes to zero, which, coupled with (3.19) gives
\begin{equation} \begin{aligned} \int _H^0 \frac {\partial }{\partial z}\left (m\frac {\partial \tilde {u}}{\partial z}\right )\mathrm{d} z= \tau _\varphi \quad \text{and} \quad \int _H^0 \frac {\partial }{\partial z}\left (m\frac {\partial \tilde {v}}{\partial z}\right )\mathrm{d} z= \tau _\theta , \end{aligned} \end{equation}
where
$\tau _\varphi \boldsymbol{e}_{\varphi } + \tau _\theta \boldsymbol{e}_{\theta }$
is the surface wind stress (cf. (3.19)). Integrating (3.16), so that it is possible to work with a baroclinic ocean without having to specify the vertical distribution of the density, and using (4.2), it follows that
\begin{align} \frac {1}{\cos \theta } \frac {\partial \varPi }{\partial \varphi }&=\tau _\varphi +2\sin \theta \, M_2, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {\partial \varPi }{\partial \theta }&=\tau _\theta -2\sin \theta \, M_1. \end{align}
Analogously, the vertically integrated continuity equation reads
\begin{equation} \frac {1}{\cos \theta }\frac {\partial M_1}{\partial \varphi } + \frac {\partial M_2}{\partial \theta } -M_2\tan \theta =0. \end{equation}
Applying
${\partial }/{\partial \theta }$
to (4.3) and
$({1}/{\cos \theta })({\partial }/{\partial \varphi })$
to (4.4), subtracting the resulting equations, and then using (4.5) and (4.3), gives
\begin{equation} 2\cos \theta \, M_2 = \frac {1}{\cos \theta }\frac {\partial \tau _\theta }{\partial \varphi }-\frac {\partial \tau _\varphi }{\partial \theta } +\tan \theta \, \tau _\varphi . \end{equation}
As discussed in Stewart (Reference Stewart2009), over large regions of the open ocean, particularly in the tropics, the wind stress is predominantly zonal, so that (in our notation)
${\partial \tau _\theta }/{\partial \varphi }\approx 0.$
Under this assumption, (4.6) reduces to
\begin{equation} 2\cos \theta M_2 \approx -\frac {\partial \tau _\varphi }{\partial \theta } +\tan \theta \,\tau _\varphi , \end{equation}
and combining this relation with (4.5) yields
\begin{equation} \frac {1}{\cos \theta }\frac {\partial M_1}{\partial \varphi } \approx \frac {1}{2\cos \theta } \left ( \frac {\partial ^2\tau _\varphi }{\partial \theta ^2} -\tan \theta \,\frac {\partial \tau _\varphi }{\partial \theta } -\sec ^2\theta \,\tau _\varphi \right )\!. \end{equation}
For comparison, the corresponding expression obtained within the tangent-plane approximation is (Sverdrup Reference Sverdrup1947; Stewart Reference Stewart2009)
\begin{equation} \frac {\partial M_x}{\partial x} \approx -\frac {1}{2\varOmega '\cos \theta } \left ( \frac {\partial \tau _x}{\partial y}\tan \theta - R'\frac {\partial ^2\tau _x}{\partial y^2} \right )\!, \end{equation}
where we have retained the notation of Sverdrup (Reference Sverdrup1947) and Stewart (Reference Stewart2009); the parameters
$\varOmega '$
and
$R'$
appear because (4.9) is written in dimensional form. The most noticeable difference between the two expressions is the presence of the term
$\sec ^2\theta \,\tau _\varphi$
in the spherical formulation. This contribution has no counterpart in the Cartesian approximation and arises solely from the exact spherical geometry retained in the present derivation.
We also note that the right-hand side of (4.6) can be recognised as the spherical curl of the wind-stress field on the unit sphere
$\mathbb{S}^2$
, and the continuity equation (4.5) implies the existence of a mass-transport streamfunction
$\varPsi$
satisfying
\begin{equation} M_1=-\frac {\partial \varPsi }{\partial \theta } \quad \text{and} \quad M_2=\frac {1}{\cos \theta }\frac {\partial \varPsi }{\partial \varphi }, \end{equation}
leading to
\begin{equation} 2\frac {\partial \varPsi }{\partial \varphi }=\text{curl}_{\mathbb{S}^2}\boldsymbol{\tau }, \quad \text{where} \quad \boldsymbol{\tau }=\tau _\varphi \boldsymbol{e}_{\varphi } + \tau _\theta \boldsymbol{e}_{\theta }. \end{equation}
The equations derived above are closely related to those obtained under the classical
$\beta$
-plane approximation (Sverdrup Reference Sverdrup1947; Stewart Reference Stewart2009). However, a crucial difference is that the present formulation fully retains the spherical geometry of the Earth. For instance, the curl operator appearing in (4.11) is not the standard Cartesian curl, but rather the curl operator on the sphere
$\mathbb{S}^2$
.
5. Ekman flows
From now on, we will restrict our attention to flows away from the equatorial region. The reason is that, for
$\theta \neq 0$
, it is convenient to perform the usual splitting of the horizontal velocity into a geostrophic and an Ekman component
\begin{equation} \begin{pmatrix} \tilde {u}\\[0.2em] \tilde {v} \end{pmatrix} = \begin{pmatrix} \tilde {u}_{{g}}\\[0.2em] \tilde {v}_{{g}} \end{pmatrix} + \begin{pmatrix} \tilde {u}_{{e}}\\[0.2em] \tilde {v}_{{e}} \end{pmatrix}\!, \end{equation}
where the geostrophic flow
$(\tilde {u}_{{g}},\tilde {v}_{{g}})$
is defined as the solution to the system
\begin{equation} \left .\begin{aligned} 2\tilde {v}_{{g}}\sin \theta &= \frac {1}{\rho \cos \theta }\frac {\partial \tilde {p}}{\partial \varphi } \\[3pt] 2\tilde {u}_{{g}}\sin \theta &= -\frac {1}{\rho }\frac {\partial \tilde {p}}{\partial \theta } \end{aligned}\;\right \} \quad \text{in} \quad \{H \lt z \lt 0\}. \end{equation}
Moreover, the equatorial region exhibits dynamical features that are absent at higher latitudes. Under the action of the trade winds, Ekman transports on either side of the equator are directed away from it, producing a strong divergence of surface waters and intense equatorial upwelling. The associated pressure and thermocline anomalies are not confined to the forcing region but are redistributed across ocean basins by equatorially trapped Kelvin and Rossby waves (see Gill Reference Gill1982 and Boyd Reference Boyd2018). As a consequence, the oceanic response cannot be described solely through a local balance between pressure gradients, Coriolis forces and wind stress. These wave-mediated adjustment processes play a central role in large-scale tropical variability, including the El Niño-Southern Oscillation. For these reasons, the equatorial band requires a distinct theoretical framework and will not be considered here.
Note that, integrating (3.16) and using (3.18), it follows that
\begin{equation} \tilde {p}(\varphi ,\theta ,z) = p_{{s}}(\varphi ,\theta ) + g\int _z^0\rho (s)\,\textrm {d}s, \end{equation}
so the geostrophic component of the velocity is uniquely determined by the surface pressure
\begin{equation} 2\tilde {v}_{{g}} \sin \theta = \frac {1}{\rho \cos \theta }\frac {\partial p_{{s}}}{\partial \varphi } \quad \text{and} \quad 2\tilde {u}_{{g}}\sin \theta = -\frac {1}{\rho }\frac {\partial p_{{s}}}{\partial \theta } \quad \text{in} \quad \{H \lt z \lt 0\}. \end{equation}
In particular, we see that, in the presence of a flat surface, the geostrophic flow vanishes when the atmospheric pressure
$p_{{s}}$
is constant. Furthermore, note that
\begin{equation} \frac {\partial \tilde {u}_{{g}}}{\partial \varphi } + \frac {\partial }{\partial \theta }(\tilde {v}_{{g}}\cos \theta ) = 0. \end{equation}
Let us introduce the complex notation
\begin{equation} W=\tilde {u}_{{e}}+\mathrm{i}\tilde {v}_{{e}}, \quad W_{{g}} = \tilde {u}_{{g}}+\mathrm{i}\tilde {v}_{{g}}, \quad \text{and} \quad W_{{w}} = u_{{a}} - \tilde {u}_{{g}} + \mathrm{i}(v_{{a}} - \tilde {v}_{{g}}). \end{equation}
The Ekman velocity
$W$
describes a boundary layer correction of the geostrophic velocity
$W_{{g}}$
due to the wind stress. We then obtain from (3.16), (3.17), (3.19) and (5.2)
\begin{equation} \left \{\begin{aligned} & \frac {\partial }{\partial z}\left (m\frac {\partial W}{\partial z}\right ) = 2\mathrm{i}\sin \theta \,\rho W - \frac {\partial }{\partial z}\left (m\frac {\partial W_{{g}}}{\partial z}\right ) & &\quad \text{in} \quad \{H \lt z \lt 0\}, \\[9pt] & \frac {\partial W}{\partial z} = -\frac {\partial W_{{g}}}{\partial z} + |W_{{w}}-\sigma W|(W_{{w}}-\sigma W) & &\quad \text{on} \quad \{z=0\}, \\[9pt] & W + W_{{g}} = 0 & &\quad \text{on} \quad \{z=H\}, \end{aligned}\right . \end{equation}
for the boundary value problem that describes the Ekman component of the flow. As we pointed out previously, if the atmospheric pressure is constant, then the geostrophic flow vanishes and (5.7) reduces to
\begin{equation} \left \{\begin{aligned} & \frac {\partial }{\partial z}\left (m\frac {\partial W}{\partial z}\right ) = 2\mathrm{i}\sin \theta \,\rho W & &\quad \text{in} \quad \{H \lt z \lt 0\}, \\[9pt] & \frac {\partial W}{\partial z} = |W_{{w}}-\sigma W|(W_{{w}}-\sigma W) & &\quad \text{on} \quad \{z=0\}, \\[9pt] & W = 0 & &\quad \text{on} \quad \{z=H\}. \end{aligned}\right . \end{equation}
Henceforth, we will focus only on the problem (5.8), leaving (5.7) to future work.
Throughout this section, to keep the notation as agile as possible, we will denote only the dependence on
$z$
explicitly unless necessary, since the other variables enter the equation only parametrically.
5.1. Some analytical results
Interestingly, the problem (5.8) is nonlinear, because even though the equation itself is linear, the upper boundary condition is nonlinear. Therefore, it is not clear a priori whether this problem is solvable and, if so, whether the solution is unique. This is the subject of the next theorem, which not only establishes the existence of a unique solution to (5.8) for given
$W_{{w}}$
,
$m$
and
$\rho$
, but also provides some information on the qualitative behaviour of the solution.
Theorem 1.
Let
$W_{{w}} \in \mathbb{C}\setminus \{0\}$
be arbitrary,
$\theta \in (-( {\pi }/{2}), {\pi }/{2})\setminus \{0\}$
, and suppose that
$\rho \in C([H,0])$
and
$m\in C^1((H,0))\cap C([H,0])$
are both positive functions. Then there exists a unique function
$W\in C^2((H,0)) \cap C([H,0])$
that solves (
5.8
). Moreover, writing
\begin{equation} W(z) = M(z)\mathrm{e}^{\mathrm{i}\vartheta (z)}, \end{equation}
where
$M(z) = |W(z)|$
, we have that, for all
$z\in (H,0)$
,
\begin{equation} M'(z) \geqslant 0 \quad \mathrm{and} \quad \mathrm{sign}(\theta )\vartheta '(z) \gt 0, \end{equation}
that is, the solution behaves as a classical Ekman spiral.
The proof can be found in Appendix A. Theorem1 remains valid if the eddy viscosity
$m$
is only piecewise continuously differentiable (but continuous everywhere), as one can see upon inspection of the proof.
It is also possible to infer a formula for the surface deflection angle, which is defined as the angle between the surface Ekman current
$W(0)$
and the difference
$W_{{w}} = W_{{a}} - W_{{g}}$
between the wind and the geostrophic current. This is the object of the next result.
Theorem 2. With
\begin{equation} \lambda = \frac {\dfrac {\partial W}{\partial z}(0)}{W(0)}, \end{equation}
let
$\chi$
be the unique positive root of the polynomial
\begin{equation} p(x) = \mu ^2x^4 + 2\mu \mathrm{Re}(\lambda )x^3 + |\lambda |^2x^2 - |\lambda W_{{w}}|^2. \end{equation}
Denote
\begin{equation} W(0) = |W(0)|\mathrm{e}^{\mathrm{i}\varTheta } \quad \mathrm{and} \quad W_{{w}} = |W_{{w}}|\mathrm{e}^{\mathrm{i}\varGamma }. \end{equation}
Then the surface deflection angle
$\omega := \varGamma - \varTheta$
is determined by the formula
\begin{equation} \tan (\omega ) = \frac {\mathrm{Im}(\lambda )}{\sigma \chi + \mathrm{Re}(\lambda )}, \end{equation}
hence, since
$\text{sign}(\theta )\mathrm{Im}(\lambda ) \gt 0$
and
$\mathrm{Re}(\lambda ) \gt 0$
, the deflection is to the right in the northern hemisphere and to the left in the southern hemisphere. In particular, we have the bounds
\begin{equation} \frac {|\mathrm{Im}(\lambda )|}{\alpha + \mathrm{Re}(\lambda )} \leq |\tan (\omega )| \leq \frac {|\mathrm{Im}(\lambda )|}{\beta +\mathrm{Re}(\lambda )}, \end{equation}
where
\begin{equation} \alpha = \frac {-\mathrm{Re}(\lambda ) + \sqrt {\mathrm{Re}(\lambda )^2 + 4\mu |\lambda W_{{w}}|}}{2} \quad \text{and} \quad \beta = \frac {-|\lambda | + \sqrt {|\lambda |^2 + 4\mu |\lambda W_{{w}}|}}{2}. \end{equation}
Moreover, although our solution is not defined for
$\theta =0$
, where the definition of the geostrophic flow breaks down, we recover the classical observation that the deflection angle vanishes as one approaches the equator (Wunsch Reference Wunsch2024); see figure 3 for a sketch.
Theorem 3.
With the notation of
Theorem 2
, for the surface deflection angle
$\omega (\theta )$
we have
\begin{equation} \lim _{\theta \to 0} \omega (\theta ) = 0. \end{equation}
Conceptual diagram showing the progressive decrease in the wind–current deflection angle toward the equator as the Coriolis force approaches zero. Not to scale.
Finally, we define the Ekman transport (denoted by
$\textrm {Ek}$
) as the integrated mass transport induced by the Ekman flow:
\begin{equation} \mathrm{Ek} = \int ^0_{H}W(z)\, \mathrm{d} z. \end{equation}
Using (5.8), we immediately see that
\begin{equation} \begin{aligned} \textrm {Ek} &= -\mathrm{i}\,\frac {m(0)}{2\sin \theta }|W_{{w}}-\sigma W(0)|(W_{{w}}-\sigma W(0)) + \underbrace {\mathrm{i}\frac {m(H)}{2\sin \theta }\frac {\partial W}{\partial z}(H)}_{\approx 0}\\ &\approx -\,\frac {\mathrm{i}}{2\sin \theta }|W_{{w}}-\sigma W(0)|(W_{{w}}-\sigma W(0)), \end{aligned} \end{equation}
since
$m(H)$
is very small and we expect
${\partial W}/{\partial z}$
to be negligible at great depths. In other words, the Ekman transport is always essentially at a right angle to the vector
$W_{{w}}-\sigma W(0)$
(or, equivalently, to the wind stress on the surface, in view of (2.9)).
5.2. Explicit solutions to the leading-order problem
This section is devoted to gaining additional insight into the qualitative behaviour of the leading-order Ekman flow by investigating some concrete examples of explicit eddy-viscosity profiles, along the lines of Constantin & Johnson (Reference Constantin and Johnson2019b
) and Zikanov et al. (Reference Zikanov, Slinn and Dhanak2003). These choices allow us to write down the solution explicitly, calculate the surface deflection angle, compare our results with some available measurements (see table 1) and visualise the flow by means of numerical plots. For simplicity, throughout this section we will assume the density to be constant:
$\rho (z) = 1$
for all
$z$
.
For clarity and ease of reference later, let us recall the leading-order problem (5.8) here, in the special case
$\rho (z) = 1$
. This problem consists of the complex-valued linear ODE
\begin{equation} \frac {\partial }{\partial z}\left (m\frac {\partial W}{\partial z}\right ) = 2\mathrm{i} W\sin \theta \quad \text{in} \quad \{H \lt z \lt 0\}, \end{equation}
with boundary conditions
\begin{equation} \frac {\partial W}{\partial z} = |W_{{w}}-\sigma W|(W_{{w}}-\sigma W) \quad \text{on} \quad \{z=0\} \end{equation}
and
\begin{equation} W = 0 \quad \text{on} \quad \{z=H\}, \end{equation}
where
$W_{{w}}\in \mathbb{C}$
corresponds to the non-dimensional wind velocity (see (5.6)) and can thus be assumed to be given.
We will consider an eddy viscosity of the form
\begin{equation} m(z) = \begin{cases} 1+\dfrac {M-1}{{z_{\star }}}\,z, &\quad z \in [{z_{\star }},0], \\[0.2em] m_1(z), &\quad z\in [z_0,{z_{\star }}], \\[0.2em] \mathfrak{m}, &\quad z \in [H,z_0], \end{cases} \end{equation}
where
$M,\mathfrak{m}\gt 0$
and the function
$m_1$
is positive, continuous and satisfies
$m_1({z_{\star }}) = M$
and
$m_1(z_0) = \mathfrak{m}$
. In other words, the eddy viscosity increases linearly to the depth
$z_{\star }$
, where it takes the value
$M$
, and is given by
$m_1$
in the interval
$[z_0,z_\star ]$
, below which it coincides with the molecular viscosity
$\mathfrak{m}$
. This increasing–decreasing profile is consistent with the so-called KPP parametrisation (see, e.g. Large et al. Reference Large, McWilliams and Doney1994; McWilliams & Huckle Reference McWilliams and Huckle2006) and, despite the absence of wave dynamics in our model, can be regarded as a simplified counterpart of the formulation proposed by Kitaigorodsky (Reference Kitaigorodsky1961) (see also the discussion by Weber Reference Weber1981). There, the eddy viscosity is taken to be proportional to the shear of the wave orbital velocity and to the square of the mixing length.
The value
$M$
can be interpreted as the maximum efficiency of turbulent momentum transport within the water column, with a larger
$M$
corresponding to more vigorous mixing and a more efficient vertical redistribution of the momentum imparted by the wind, while the depth
$z_{\star }$
characterises the depth at which the wind-generated turbulence is at its strongest. We will return to this important fact in § 6.
Recalling that we are setting
$m(0)=1$
, it can be assumed that
$M\approx 10$
(see McWilliams & Huckle Reference McWilliams and Huckle2006 and McWilliams et al. Reference McWilliams, Huckle, Liang and Sullivan2012). In some works (see, e.g. Weber Reference Weber1981; Zikanov et al. Reference Zikanov, Slinn and Dhanak2003), the value of
$M$
appears to be greater than
$10$
, even if its value is not easily quantifiable. In the latter case, we will set
$M=50$
, in order to analyse how a higher value of
$M$
affects the Ekman flow. Moreover, we will set
$\mathfrak{m} = 0.01$
for the molecular viscosity. Finally, according to Zikanov et al. (Reference Zikanov, Slinn and Dhanak2003), McWilliams & Huckle (Reference McWilliams and Huckle2006) and McWilliams et al. (Reference McWilliams, Huckle, Liang and Sullivan2012), one can estimate
$z_{\star }\approx z_0/4$
,
$z_{\star }\approx z_0/5$
and
$z_{\star }\approx z_0/6$
, respectively, with
$z_0\approx 5$
(corresponding to
$60$
m) according to McWilliams & Huckle (Reference McWilliams and Huckle2006) and
$z_0\approx 18$
(corresponding to
$200$
m) according to McWilliams et al. (Reference McWilliams, Huckle, Liang and Sullivan2012). For particular choices of
$m_1$
, the equations can be solved explicitly; see Appendix B for a brief discussion of the general solving method in this case.
In the sequel, we will consider some of such situations and compare the results with available measurements. First, we consider the easiest possible case, that is, that of a constant
$m_1$
. To avoid dealing with discontinuous eddy viscosity, we take
$z_0=H$
in this case; while this example is physically of dubious relevance, it is nevertheless interesting to make a comparison with a model similar to that of Ekman (Reference Ekman1905). Next, we investigate the cases of linearly decaying eddy viscosity, which can be identified as an approximation of the viscosity profiles in Zikanov et al. (Reference Zikanov, Slinn and Dhanak2003), McWilliams & Huckle (Reference McWilliams and Huckle2006) and McWilliams et al. (Reference McWilliams, Huckle, Liang and Sullivan2012). Finally, we turn our attention to the exponentially decaying eddy viscosity proposed by Wenegrat & McPhaden (Reference Wenegrat and McPhaden2016) to match the measurements by Peters et al. (Reference Peters, Gregg and Toole1988) and Dillon et al. (Reference Dillon, Moum, Chereskin and Caldwell1989), where the solution can be written in terms of modified Bessel functions. For details on the computations involved, see Appendix B.
5.2.1.
Constant
$m_1$
The first case we consider is that of a constant
$m_1$
, where we take only two layers (that is,
$z_0=H$
)
\begin{equation} m_1(z) = M, \quad z\in [H,{z_{\star }}]. \end{equation}
Although the assumption of constant eddy viscosity may seem unrealistic, it has been observed that it may provide a good fit to certain wind-driven drift currents (Polton et al. Reference Polton, Lenn, Elipot, Chereskin and Sprintall2013). In this case, the solution is explicit (see Appendix B1) and is depicted in figure 4 (see also figure 11 in Appendix D). Figure 5 presents the computed deflection angles. All values remain below
$45^\circ$
, owing to the linear portion of the eddy-viscosity profile adopted in the model.
The velocity field
$W = \tilde {u}_{{e}} + \mathrm{i}\tilde {v}_{{e}}$
in the case of constant
$m_1$
for
$|W_{{w}}|=1$
,
$\sigma = 0.014$
(corresponding to
$U_{ {a}}' = 10\,\mathrm{m\,s^{-1}}$
),
$\theta =45^\circ$
and
$H=-80$
. In all cases, the plots are displayed only up to
$z=-30$
: (a)
$M = 10$
,
${z_{\star }} = -1$
, (b)
$M = 50$
,
${z_{\star }} = -1$
.
The surface deflection angle in the case of constant
$m_1$
for
$|W_{{w}}|=1$
,
$\sigma = 0.014$
and
$H=-80$
: (a)
$M = 10$
, (b)
$M = 50$
.
5.2.2.
Linearly decaying
$m_1$
The second case we study is of an eddy viscosity that decays linearly in the lower layer to a small value, denoted by
$\mathfrak{m}$
\begin{equation} m_1(z) = \frac {{z_{\star }}\mathfrak{m}-z_0 M}{{z_{\star }}-z_0} + \frac {M-\mathfrak{m}}{{z_{\star }}-z_0}\,z, \quad z\in [z_0,{z_{\star }}]. \end{equation}
Since the concept of eddy viscosity is inherently associated with the presence of turbulence, which in the present case is generated by wind forcing, it is reasonable to assume that below the Ekman layer, where wind-induced effects are negligible, turbulence should be absent and the viscosity for
$z\lt z_0$
should reduce the molecular viscosity, which is typically several orders of magnitude smaller than the maximum value of the eddy viscosity; therefore, we set
$\mathfrak{m}=0.01$
. See figures 6 and 7 for plots of the solution and the deflection angle for different values of the parameters (as well as figure 12 in Appendix D).
The velocity field
$W = \tilde {u}_{{e}} + \mathrm{i}\tilde {v}_{{e}}$
in the case of linearly decaying
$m_1$
for
$|W_{{w}}|=1$
,
$\sigma = 0.014$
,
$H=-80$
and
$\mathfrak{m}=0.01$
. In all cases, the plots are displayed only up to
$z=-21$
: (a)
$M=10$
,
$z_0=-5$
,
$z_{\star }=-1$
, (b)
$M=10$
,
$z_0=-20$
,
$z_{\star }=-4$
(c)
$M=50$
,
$z_0=-5$
,
$z_{\star }=-1$
, (d)
$M=50$
,
$z_0=-20$
,
$z_{\star }=-4$
.
Figure 6. Long description
Panel A: Two line graphs show the velocity components ue and ve as functions of depth z. The left graph plots ue with z ranging from -21 to 0 and ue from -0.2 to 0.6. The right graph plots ve with z ranging from -21 to 0 and ve from -0.4 to 0. Panel B: Two line graphs show the velocity components ue and ve as functions of depth z. The left graph plots ue with z ranging from -21 to 0 and ue from -0.2 to 0.6. The right graph plots ve with z ranging from -21 to 0 and ve from -0.4 to 0. Panel C: Two line graphs show the velocity components ue and ve as functions of depth z. The left graph plots ue with z ranging from -21 to 0 and ue from -0.2 to 0.6. The right graph plots ve with z ranging from -21 to 0 and ve from -0.4 to 0. Panel D: Two line graphs show the velocity components ue and ve as functions of depth z. The left graph plots ue with z ranging from -21 to 0 and ue from -0.2 to 0.6. The right graph plots ve with z ranging from -21 to 0 and ve from -0.4 to 0.
The surface deflection angle in the case of linearly decaying
$m_1$
for
$|W_{{w}}|=1$
,
$\sigma = 0.014$
,
$H=-80$
and
$\mathfrak{m}=0.01$
: (a)
$M=10$
,
$z_0=-5$
, (b)
$M=10$
,
$z_0=-20$
(c)
$M=50$
,
$z_0=-5$
(d)
$M=50$
,
$z_0=-20$
.
One remark is in order. The velocity profiles shown in figures 6(a) and 6(c), although differentiable everywhere, exhibit noticeably more irregular behaviour than those in figure 6(b). This can be attributed to the eddy-viscosity profile and, in particular, to its degree of smoothness. For the piecewise linear eddy viscosity considered here, the discontinuity in the derivative of
$m_1$
at
$z_0$
, namely
${M-\mathfrak{m}}/{z_\star -z_0}$
, is substantially larger in the cases shown in figures 6(a) and 6(c) than in figure 6(b). The latter therefore corresponds to the smoothest velocity profiles among those displayed in figure 6. This effect becomes even more apparent when comparing figure 6 with figure 8, which shows the velocity profiles obtained with an exponentially decaying eddy viscosity. In that case, the jump in the derivative of
$m_1$
at
$z_0$
is given by
$\mathfrak{m}( {\ln (\mathfrak{m}/M )}/{z_0-z_\star })$
. For example, with
$M=50$
,
$z_0=-5$
,
$z_\star =-1$
and
$\mathfrak{m}=0.01$
, the jump is approximately
$12.5$
for the piecewise linear eddy viscosity (cf. figure 6
c), whereas it is only approximately
$0.02$
for the exponentially decaying profile. This substantial reduction in the derivative discontinuity explains the markedly smoother velocity profiles observed in figure 8.
5.2.3. Exponentially decaying eddy viscosity
The velocity field
$W = \tilde {u}_{{e}} + \mathrm{i}\tilde {v}_{{e}}$
in the case of exponentially decaying
$m_1$
with
$|W_{{w}}|=1$
,
$\sigma = 0.014$
. Here,
$z_0=-5$
,
$H=-80$
and
$\mathfrak{m}=0.01$
. In all cases, the plots are displayed only up to
$z=-12$
: (a)
$M=10$
,
$z_0=-5$
,
$z_\star =-1$
, (b)
$M=10$
,
$z_0=-20$
,
$z_\star =-4$
, (c)
$M=50$
,
$z_0=-5$
,
$z_\star =-1$
, (d)
$M=50$
,
$z_0=-20$
,
$z_\star =-4$
.
Figure 8. Long description
The image contains four sets of line graphs, each set consisting of two graphs labeled (a), (b), (c), and (d). Each graph shows the relationship between depth (z) and two velocity components (u_e and v_e). Panel A: For M = 10, z0 = -5, z* = -1, the left graph shows u_e versus z, and the right graph shows v_e versus z. The u_e graph has a range from -0.2 to 0.6, while the v_e graph ranges from -0.4 to 0. Panel B: For M = 10, z0 = -20, z* = -4, the left graph shows u_e versus z, and the right graph shows v_e versus z. The u_e graph ranges from -0.2 to 0.6, and the v_e graph ranges from -0.4 to 0. Panel C: For M = 50, z0 = -5, z* = -1, the left graph shows u_e versus z, and the right graph shows v_e versus z. The u_e graph ranges from -0.2 to 0.6, while the v_e graph ranges from -0.4 to 0. Panel D: For M = 50, z0 = -20, z* = -4, the left graph shows u_e versus z, and the right graph shows v_e versus z. The u_e graph ranges from -0.2 to 0.6, and the v_e graph ranges from -0.4 to 0. Each graph displays how the velocity components change with depth for different parameter values.
The final case we consider is that proposed as a model dictated by some turbulence measurements in non-equatorial regions (Cronin & Kessler Reference Cronin and Kessler2009; Wenegrat & McPhaden Reference Wenegrat and McPhaden2016; Sentchev et al. Reference Sentchev, Yaremchuk, Bourras, Pairaud and Fraunié2023) of an exponentially decreasing eddy viscosity of the form
\begin{equation} m_1(z)=M\mathrm{e}^{q(z-{z_{\star }})}, \quad z_0 \lt z \lt 0, \end{equation}
where
\begin{equation} q=\frac {\ln \left (\dfrac {\mathfrak{m}}{M}\right )}{z_0-{z_{\star }}}, \end{equation}
which ensures that
$m(z_0)=\mathfrak{m}$
. Analogously to the previous case, we set
$\mathfrak{m}=0.01$
. The solution and the deflection angle are shown in figures 8 and 9 for different values of the parameters; see also figure 13 in Appendix D. Observe that the velocity profiles in figure 8 are essentially a smoother version of those in figure 6.
The surface deflection angle in the case of exponentially decaying
$m_1$
for
$|W_{{w}}|=1$
,
$\sigma = 0.014$
,
$H=-80$
and
$\mathfrak{m}=0.01$
: (a)
$M=10$
,
$z_0=-5$
, (b)
$M=10$
,
$z_0=-20$
, (c)
$M=50$
,
$z_0=-5$
, (d)
$M=50$
,
$z_0=-20$
.
6. Discussion
We conclude with a summary of what this work has achieved and an overview of the possible directions for future research that it opens up.
Starting from the governing equations for a viscous fluid in spherical coordinates, we derived an asymptotic model for the analysis of large-scale wind-drift currents. Unlike classical approximations such as the
$f$
-plane or
$\beta$
-plane, our approach retains the full spherical geometry, in the spirit of recent works by Constantin & Johnson (Reference Constantin and Johnson2018, Reference Constantin and Johnson2019a
,Reference Constantin and Johnson
b
), and provides a consistent framework for investigating wind-driven flows on large spatial scales, with direct relevance to ocean circulation (Wunsch & Ferrari Reference Wunsch and Ferrari2004). The derivation was based on asymptotic expansions with respect to the two small parameters
\begin{equation} \varepsilon =\frac {D'}{R'} \quad \text{and} \quad \mathscr{R}=\frac {U'}{\varOmega 'R'}, \end{equation}
which characterise the problem and whose smallness reflects the underlying physical regime. Importantly, both parameters emerge naturally from the non-dimensionalisation of the governing equations and do not require any additional assumptions. The scales
$R'$
and
$\varOmega '$
arise directly from the geometry and rotation of the Earth. The characteristic depth scale
$D'$
reflects the assumption that the relevant flow occurs near the ocean surface; this leads to the standard scaling in which the vertical velocity is of order
$\varepsilon$
smaller than the horizontal velocity components, implying a weak coupling between vertical and horizontal motions (see, for instance, Constantin & Johnson Reference Constantin and Johnson2019a
,
Reference Constantin and Johnsonb
and Puntini Reference Puntini2025, Reference Puntini2026). The Rossby number
$\mathscr{R}$
measures the relative importance of inertial and rotational effects; its smallness identifies a regime in which the Earth’s rotation dominates over advection.
The asymptotic expansion provided here highlights the fact that, for such flows, the linear terms are the ones determining – at leading order – the dynamics of the wind-driven currents. This theoretical result is related to the fact that, as discussed by Ferrari & Wunsch (Reference Ferrari and Wunsch2010) and Wunsch (Reference Wunsch2011), the large-scale motion is deeply connected with the distribution of the potential energy of the ocean, while the predominant contribution to the kinetic energy comes from transport and mixing at much smaller spatial scales; moreover, the time scale for injection/dissipation of kinetic energy (stemming from the nonlinear terms in the equations of motion) is much smaller than the one in (2.16), that is, the natural choice
$T' = R'/U'$
. Moreover, recent studies such as those by Wunsch (Reference Wunsch2023, Reference Wunsch2024), performing time averages over many years, argue that there are some simple descriptive, but quantitative, pattern properties of the ocean circulation that are near-globally applicable. Our theoretical analysis can be, in a sense, viewed as a counterpart of these data-driven analyses.
Only two assumptions were made a priori: turbulent stresses are parametrised through an eddy-viscosity closure, and the eddy viscosity and the density are assumed to depend solely on depth. As we pointed out at the beginning of § 2.3, the eddy-viscosity formulation comes with a caveat. In the present work, we assume that, on the surface
$\{z=0\}$
, the eddy viscosity attains a finite value that is subsequently used as the reference scale for the non-dimensionalisation. However, from a physical point of view, the notion of a local eddy viscosity at the ocean surface is not entirely unproblematic. Near the air–sea interface, turbulence is strongly influenced by wave breaking and other air–sea interaction processes, which can induce non-local momentum transport. Consequently, the turbulent stress cannot, in general, be related to the local mean velocity shear through a simple diffusive (eddy viscosity) closure (see, e.g. Craig & Banner Reference Craig and Banner1994). The surface value of the eddy viscosity should therefore be regarded primarily as a convenient scaling parameter rather than as a quantity with a direct physical interpretation.
Later on, we made a further simplifying assumption, namely that the free surface should be flat. Although, as we pointed out earlier, this assumption can be justified in view of the weak coupling between vertical and horizontal velocities, it is nonetheless an obvious shortcoming, since it neglects physical processes that are known to play a central role in the dynamics of the upper ocean. In particular, the drift in the upper layer is not driven solely by mean ocean currents, but is also influenced by surface waves (see, for example, Rascle & Ardhuin Reference Rascle and Ardhuin2009; Röhrs et al. Reference Röhrs, Christensen, Hole, Broström, Drivdal and Sundby2012), and, in the upper few metres of the ocean, the wave-induced Stokes drift contributes significantly to the overall transport and structure of the flow (Stokes Reference Stokes1847; Monismith & Fong Reference Monismith and Fong2004; Röhrs et al. Reference Röhrs, Christensen, Vikebø, Sundby, Saetra and Broström2014). This influence is further enhanced by the dissipation of wave momentum (Carniel et al. Reference Carniel, Warner, Chiggiato and Sclavo2009; Christensen & Terrile Reference Christensen and Terrile2009), which couples with the Earth’s rotation through the Coriolis–Stokes force. The resulting wave–current interactions can substantially modify both the magnitude and direction of near-surface currents (Ursell Reference Ursell1950; Hasselmann Reference Hasselmann1970; Lewis & Belcher Reference Lewis and Belcher2004; Polton et al. Reference Polton, Lewis and Belcher2005). However, in principle, our model can also account for wave motion: the kinematic boundary (3.5) provides an equation for the free surface, which also enters the equations directly through the dynamic boundary conditions (2.23). In particular, denoting by
$z = \tilde {h}$
the free surface at leading order, (5.3) would read
\begin{equation} \tilde {p} = p_{{s}} + g\int _z^{\tilde {h}}\rho (s)\,\textrm {d}s, \end{equation}
and thus the geostrophic flow would be given by
\begin{equation} \left .\begin{aligned} 2\tilde {v}_{{g}}\sin \theta &= \frac {1}{\rho \cos \theta }\left (\frac {\partial p_{{s}}}{\partial \varphi } + g\rho (\tilde {h})\frac {\partial \tilde {h}}{\partial \varphi }\right ) \\ 2\tilde {u}_{{g}}\sin \theta &= -\frac {1}{\rho }\frac {\partial p_{{s}}}{\partial \theta } - g\frac {\rho (\tilde {h})}{\rho }\frac {\partial \tilde {h}}{\partial \theta } \end{aligned}\;\right \} \quad \text{in} \quad \{H \lt z \lt \tilde {h}\}. \end{equation}
Clearly, this significantly complicates the analysis, even in the case of constant atmospheric pressure; this is a problem that we hope to tackle in future work.
Even with the simplification of a flat surface, our model has interesting features. One of the key aspects is that it does not rely on tangent-plane approximations; instead, it fully retains the effects of the sphericity of the Earth. To exemplify the potential of this formulation, in § 4 we discussed the simplest model of ocean circulation, that is, Sverdrup’s model (Sverdrup Reference Sverdrup1947). Already in this simplified context, the formulation in global spherical coordinates gives rise to additional terms in the formulae due to the variations of
$\theta$
. We believe that this approach could be applied successfully also to less simple models and provide a more accurate description of large-scale flows than is possible within the framework of tangent-plane theory; this is something that we plan to pursue in later work.
In § 5, we focused on Ekman flows. There, we split the velocity field into a geostrophic and an Ekman component (the latter being essentially a boundary layer correction that takes into account the effects of viscosity close to the flat surface). It should be noted that this splitting requires the latitude
$\theta$
to be non-zero (that is, that we are away from the equator), because the definition (5.2) of the geostrophic flow requires dividing by
$\sin \theta$
. It is also to be expected that for small values of
$\theta$
(say, for
$|\theta |$
below
$5^\circ$
) this approach will not be an accurate description of the real flow, because close to the equator the Coriolis force becomes small and the linear Coriolis terms do not dominate the nonlinear advective terms in the Navier–Stokes equations (contrarily to what happens at higher latitudes), so that a different scaling as in Stommel’s theory of wind-driven drift currents (Stommel Reference Stommel1960; Webb Reference Webb2018) would be required to capture the emergence of nonlinear effects at leading order; see also the remarks at the beginning of § 5.
Keeping these considerations in mind and assuming for simplicity that the atmospheric pressure is constant (which, together with the assumption of a flat surface, implies that the underlying geostrophic flow vanishes – cf. (6.3)), we proved analytically that the leading-order problem possesses a unique solution that behaves qualitatively like an Ekman spiral, and we were able to characterise the surface deflection angle and the intensity of the surface current, showing in particular that the deflection angle goes to zero as the latitude
$\theta$
goes to zero. Then, to obtain some explicit values to compare with earlier solutions and with experimental data available in the literature, we assumed the density to be constant (to allow us to look for explicit solutions) and adopted three different explicit parametrisations for the eddy viscosity. In all cases, the eddy viscosity is taken to decrease with depth from the surface until a certain depth
$z_{\star }$
, below which we have three profiles: constant, and linearly and exponentially decreasing to a depth
$z_0$
, below which the eddy viscosity becomes the constant (and very small) molecular viscosity. Choosing reasonable values for the maximum of the eddy viscosity and for the depths
$z_{\star }$
and
$z_0$
, we plotted various profiles of the Ekman spiral and, more importantly, calculated the deflection angles. These can be compared with measurements: in practice, the deflection angle can be observed to be both greater (Chereskin Reference Chereskin1995; Niiler & Paduan Reference Niiler and Paduan1995; Yoshikawa & Masuda Reference Yoshikawa and Masuda2009; Kim et al. Reference Kim, Cornuelle and Terrill2010; Yoshikawa et al. Reference Yoshikawa, Masuda, Marubayashi and Ishibashi2010; Röhrs & Christensen Reference Röhrs and Christensen2015) and lower (Poulain et al. Reference Poulain, Gerin, Mauri and Pennel2009; Yoshikawa & Masuda Reference Yoshikawa and Masuda2009; Yoshikawa et al. Reference Yoshikawa, Masuda, Marubayashi and Ishibashi2010; Sentchev et al. Reference Sentchev, Forget and Fraunié2017; Berta et al. Reference Berta2018) than 45
$^\circ$
; see the values in table 1. We recover both regimes, obtaining highly realistic values, as can be seen comparing figures 4, 6 and 8 with table 1.
Overview of reported measurements of the surface deflection angle of wind-drift currents.
Table 1. Long description
Table with multiple columns and rows summarizing measurements of the surface deflection angle of wind-drift currents. Columns include references to studies, specific measurements, and deflection angles. Each row represents a different study or measurement.
Away from equatorial regions (where our model ceases to be valid), our results from § 5.2 seem to point to the following conclusions. For fixed parameters of the eddy viscosity, in most cases the deflection angle decreases with latitude; moreover, as the maximum
$M$
of the eddy viscosity increases, the deflection angle increases, and as the depth
$z_{\star }$
of the maximum
$M$
decreases, the deflection angle increases as well; and finally, the faster the eddy viscosity decays in the interval
$[z_0,z_{\star }]$
, the higher the deflection angle.
As we mentioned in § 5.2, from a physical point of view,
$M$
represents the maximum efficiency of turbulent momentum transport within the water column. A larger
$M$
indicates more vigorous mixing, meaning that the momentum imparted by the wind is redistributed vertically more effectively. The depth
$z_{\star }$
at which this maximum occurs characterises how deeply the wind-generated turbulence penetrates before reaching its greatest intensity. For a fixed
$M$
, a larger
$z_{\star }$
implies that turbulent motions remain energetic to greater depths, allowing the influence of the surface forcing to extend further into the ocean interior. On the other hand, a smaller
$z_{\star }$
is often associated with a stronger stratification, which suppresses vertical motions and confines the turbulence closer to the surface. In this sense,
$z_{\star }$
can also be viewed as an indicator of the balance between wind-driven turbulent production and stratification-induced damping. In particular, we infer that larger deflection angles, such as those observed by Chereskin (Reference Chereskin1995) and Röhrs & Christensen (Reference Röhrs and Christensen2015), are attained for high values of
$M$
, which corresponds to a more efficient mixing due to the wind, coupled with a strong stratification (i.e. a small
$z_{\star }$
). In contrast, small deflection angles, such as those measured by Poulain et al. (Reference Poulain, Gerin, Mauri and Pennel2009), Sentchev et al. (Reference Sentchev, Forget and Fraunié2017) and Berta et al. (Reference Berta2018), arise for smaller
$M$
and weaker stratification (i.e. larger
$z_{\star }$
). Additionally, the velocity profiles provided by Sentchev et al. (Reference Sentchev, Forget and Fraunié2017) can be seen to be qualitatively similar to those of figures 6(b) and 8(b), which indeed correspond to small
$M$
and large
$z_\star$
(
$M=10$
,
$z_\star =4$
).
The case of exponentially decaying
$m_1$
with
$|W_{{w}}|=1$
,
$\sigma = 0.014$
. Here,
$z_0=-5$
,
$H=-80$
and
$\mathfrak{m}=0.01$
. In all cases, the plots are displayed only up to
$z=-5$
: (a)
$M=10$
,
$z_\star =-0.1$
, (b)
$M=10$
,
$z_\star =-1$
, (c)
$M=50$
,
$z_\star =-0.1$
, (d)
$M=50$
,
$z_\star =-1$
.
Additionally, measurements show that the speed of the Ekman current decays with depth more rapidly than its directional rotation, in contrast to the theoretical predictions from Ekman’s classical theory (Wijffels et al. Reference Wijffels, Firing and Bryden1994). This effect is known as Ekman spiral flattening or compression (Schudlich & Price Reference Schudlich and Price1998). Although it can be easily observed in our figures 4, 6 and 8, it is interesting to investigate how the Ekman flattening is affected by the parameters
$M$
and
$z_\star$
. Figure 10 shows the velocity profiles for different two different values of
$M$
and two different values
$z_\star$
, for an exponentially decaying
$m_1.$
The plots indicate that the flattening is enhanced as both
$M$
and
$z_\star$
decrease, corresponding to a reduced mixing and a higher stratification (or, in general, a turbulent motion confined close to the surface).
We conclude with one last remark. In the present paper, we have dealt with the currents generated by the wind over a free surface, namely considering air–sea interactions via the boundary condition (2.11). However, a similar analysis could be carried out for the bottom boundary layer or the atmospheric boundary layer. In this case, the ODE (5.7) remains unchanged, whereas the boundary conditions must be modified (see, e.g. Veronis Reference Veronis1980 or Stefanescu Reference Stefanescu2024). Specifically, for atmospheric flows over a rigid surface at
$\{z=0\}$
, one obtains
\begin{equation} \left \{\begin{aligned} & \frac {\partial }{\partial z}\left (m\frac {\partial W}{\partial z}\right ) = 2\mathrm{i}\rho W\sin \theta - \frac {\partial }{\partial z}\left (m\frac {\partial W_{{g}}}{\partial z}\right ) \quad \text{in} \quad \{z\gt 0\}, \\[0.4em] & W = -W_{{g}} \quad \text{on} \quad \{z=0\}, \\[0.4em] & \lim _{z\rightarrow +\infty } W(z)=0. \end{aligned}\right . \end{equation}
In the simplest setting, where the geostrophic current
$W_{{g}}$
is constant and both the eddy viscosity and the density are uniform, i.e.
$m(z)=\rho (z)=1$
, the solution to (6.4) is
\begin{equation} W(z)=-W_{{g}}\mathrm{e}^{-\sqrt {|\textrm{sin} \theta |}\,z}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\sqrt {|\textrm{sin} \theta |}\,z}. \end{equation}
However, a comprehensive discussion of boundary layer flows near a rigid surface lies beyond the scope of this work.
Acknowledgements
The authors are grateful to the anonymous reviewers for their valuable and insightful comments and suggestions that greatly helped to improve the quality of the paper. C.P. would like to thank K. Emanuel and R. Rotunno for the interesting discussion on this topic during the workshop ‘Vortices and vorticity in geophysical flows’ within the scientific trimester ‘Mathematical Developments in Geophysical Fluid Dynamics’, held at IHP in Paris.
Funding
C.P. was partially supported by the Austrian Science Fund (FWF) (grant number Z 387-N, grant doi: https://doi.org/10.55776/Z387).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
All data used in this work are properly cited throughout the manuscript.
Author contributions
C.P. proposed the topic, developed the governing equations in §§ 2 and 3, carried out the theory of § 4 as well as the theory and numerics of § 5.2 and contributed to writing, alongside minor contributions to the theory of § 5.1. L.R. developed the governing equations in §§ 2 and 3 and the theory of § 5.1 and contributed to the theory and numerics of § 5.2 and to writing. E.S. contributed to Appendix A. The revision was carried out by C.P. and L.R.
Appendix A. Proofs of the results from § 5.1
We provide here the proofs of the analytical results we formulated in § 5.1.
A.1. Proof of Theorem 1
For simplicity, we denote a derivative with respect to
$z$
by a prime. Let
$Y\in C^2((H,0))\cap C([H,0])$
be the solution to the linear initial value problem
\begin{equation} \left \{\begin{aligned} & (m(z)Y'(z))' - 2\mathrm{i}\sin \theta \,\rho (z)Y(z) = 0, \quad H \lt z \lt 0, \\[0.2em] & Y(H) = 0, \\[0.2em] & Y'(H) = 1; \end{aligned}\right . \end{equation}
the existence of a unique global solution to this problem follows from standard theory for second-order linear systems of ODEs (see, e.g. Teschl Reference Teschl2012). Multiplying the first equation of (A1) with the complex conjugate
$\overline {Y}$
of
$Y$
, integrating from
$H$
to
$z$
, performing an integration by parts and taking the third equation in (A1) into account, we obtain
\begin{equation} m(z)Y'(z)\overline {Y(z)} = \int _{H}^z m(s)|Y'(s)|^2\,\textrm {d}s + 2\mathrm{i}\sin \theta \int _{H}^z\rho (s)|Y(s)|^2\,\textrm {d}s. \end{equation}
Writing
\begin{equation} Y(z) = r(z)\mathrm{e}^{\textrm{i}\phi (z)}, \end{equation}
plugging into (A2) and taking the real and imaginary part, we see that, on the one hand,
\begin{equation} m(z)r'(z)r(z) = \int _{H}^z m(s)|Y'(s)|^2\,\textrm {d}s \gt 0, \end{equation}
hence, since
$m(z) \gt 0$
and
$r(z) \ge 0$
, we have
\begin{equation} r(z) \gt 0 \quad \text{and} \quad r'(z) \gt 0 \quad \text{for all $z\gt H$}, \end{equation}
and, on the other hand,
\begin{equation} m(z)\phi '(z)r(z)^2 = 2\sin \theta \int _{H}^z\rho (s)|Y(s)|^2\,\textrm {d}s \begin{cases} \gt 0 & \text{if $\theta \in \left (0,\dfrac {\pi }{2}\right )$,} \\[0.7em] \lt 0 & \text{if $\theta \in \left (-\dfrac {\pi }{2},0\right )$,} \end{cases} \end{equation}
thus
\begin{equation} \text{sign}(\theta )\phi '(z) \gt 0 \quad \text{for all $z \gt H$}. \end{equation}
In particular, we have
\begin{equation} Y(z) \neq 0 \quad \text{for all $z \gt H$}. \end{equation}
Now, note that a solution
$W$
to (5.8) must necessarily be of the form
\begin{equation} W = cY \end{equation}
for some
$c\in \mathbb{C}$
; therefore, we only have to show that we can find a unique
$c$
such that the nonlinear upper boundary condition is met. To this end, let us simplify the notation by setting
\begin{equation} a := W(0) \quad \text{and} \quad d := W_{{w}} -\sigma a, \end{equation}
hence
$c = {a}/{Y(0)}$
(as we saw,
$Y(0) \neq 0$
) and
\begin{equation} W'(0) = \lambda a, \end{equation}
where
\begin{equation} \lambda := \frac {Y'(0)}{Y(0)} = \frac {r'(0)}{r(0)} + \mathrm{i}\phi '(0). \end{equation}
Note that
$\mathrm{Re}(\lambda ) \gt 0$
and
$\text{sign}(\theta )\mathrm{Im}(\lambda ) \gt 0$
, in view of (A5) and (A7).
In order for the second equation in (5.8) to be satisfied, we must have
\begin{equation} \frac {\lambda }{\mu }(W_{{w}}-d) = |d|d, \end{equation}
which, writing
\begin{equation} d = \chi \mathrm{e}^{\textrm{i}\delta } \end{equation}
for
$\chi = |d|$
, becomes
\begin{equation} \lambda W_{{w}} = \chi \mathrm{e}^{\textrm{i}\delta }(\lambda + \sigma \chi ). \end{equation}
Multiplying both sides of this equation by the respective complex conjugate and rearranging terms, we see that
$\chi$
is determined through the positive solutions of the equation
\begin{equation} p(\chi ) = 0, \end{equation}
where
$p$
is the real fourth-order polynomial
\begin{equation} p(x) = \mu ^2x^4 + 2\mu \mathrm{Re}(\lambda )x^3 + |\lambda |^2x^2 - |\lambda W_{{w}}|^2. \end{equation}
Note that
$p(0) \lt 0$
,
$p(x) \to \infty$
as
$x\to \infty$
and
$p$
is strictly increasing on
$(0,\infty )$
, because
\begin{equation} p'(x) = 2x\bigl (2\mu ^2 x^2 + 3\mu \mathrm{Re}(\lambda )x + |\lambda |^2\bigr ) \gt 0 \quad \text{for $x\gt 0$}; \end{equation}
therefore, there is a unique positive solution
$\chi$
of (A16), from which, in view of (A15), we recover the constant
$c$
through
\begin{equation} c = \frac {a}{Y(0)} = \frac {W_{{w}}-d}{Y(0)} = \frac {\chi }{\sigma \chi + \lambda }\boldsymbol{\cdot }\frac {W_{{w}}}{Y(0)} \end{equation}
(note that
$\sigma \chi + \lambda \neq 0$
, since
$\chi \gt 0$
and
$\lambda \in \mathbb{C}\setminus \mathbb{R}$
). Finally, since
$W$
is a multiple of
$Y$
, the claimed properties (5.10) of the Ekman spiral follow immediately from (A5) and (A7).
A.2. Proof of Theorem 2
By (A19), we have
\begin{equation} |W(0)|\mathrm{e}^{\textrm{i}\varTheta } = W(0) = \frac {\chi }{\sigma \chi + \lambda }W_{{w}} = \frac {1}{\mu }\left (1+\frac {\lambda }{\sigma \chi }\right )^{-1}|W_{ { w}}|\mathrm{e}^{\textrm{i}\varGamma }, \end{equation}
thus
\begin{equation} \left |\frac {W_{{w}}}{W(0)}\right |\mathrm{e}^{{i}(\varGamma - \varTheta )} = \sigma + \frac {\lambda }{{\chi }}, \end{equation}
whence (5.14) follows immediately. To prove the bounds (5.15), note that, with
$p$
as in (5.12), since
$\mathrm{Re}(\lambda ) \leq |\lambda |$
we have
\begin{equation} q(x) \leq p(x) \leq r(x), \quad x\gt 0, \end{equation}
where
\begin{align} r(x) = |\lambda |^2 x^2\left (\frac {\mu }{|\lambda |}x + 1\right )^2 - |\lambda W_{{w}}|^2 \quad \text{and} \quad q(x) = \mathrm{Re}(\lambda )^2 x^2\left (\frac {\mu }{\textrm { Re}(\lambda )}x + 1\right )^2 - |\lambda W_{{w}}|^2. \end{align}
Therefore, the unique positive root
$\chi$
of
$p$
must be smaller than the unique positive root of
$q$
and larger than the unique positive root of
$r$
, which are given by
$\alpha /{\sigma }$
and
$\beta /{\sigma }$
, respectively, and the bounds (5.15) follow.
A.3. Proof of Theorem 3
Let
$Y_\theta$
be the solution to (A1) (where we explicitly denote the
$\theta$
-dependence via the index
$\theta$
) and
$Y_0$
the solution to
\begin{equation} \left \{\begin{aligned} & (m(z)Y_0'(z))' = 0, \quad H \lt z \lt 0, \\[0.2em] & Y_0(H) = 0, \\[0.2em] & Y_0'(H) = 1, \end{aligned}\right . \end{equation}
which is easily found to be
\begin{equation} Y_0(z) = \int _H^z\frac {m(H)}{m(s)}\,\textrm {d}s. \end{equation}
Then, the function
$\hat {Y}_\theta = Y_\theta -Y_0$
solves the problem
\begin{equation} \left \{\begin{aligned} & (m(z)\hat {Y}_\theta '(z))' = 2\mathrm{i}\sin \theta \,\rho (z)\hat {Y}_\theta (z) + 2\mathrm{i}\sin \theta \,\rho (z)Y_0(z), \quad H \lt z \lt 0, \\[0.2em] & \hat {Y}_\theta (H) = 0, \\[0.2em] & \hat {Y}_\theta '(H) = 0. \end{aligned}\right . \end{equation}
Integrating the first equation yields
\begin{equation} \hat {Y}_\theta '(z) = \frac {2\mathrm{i}\sin \theta }{m(z)}\int _H^z\rho (s)\hat {Y}_\theta (s)\,\textrm {d}s + \frac {2\mathrm{i}\sin \theta }{m(z)}\int _H^z\rho (s)Y_0(s)\,\textrm {d}s; \end{equation}
integrating once more, changing variables and taking absolute values, we obtain
\begin{equation} |\hat {Y}_\theta (z)| \leq 2|\sin \theta |\int _H^z\left (\int _t^0\frac {\mathrm{d} s}{m(s)}\right )\rho (t)|\hat {Y}_\theta (t)|\,\textrm {d}t + 2|\sin \theta |\int _H^z\left (\int _t^z\frac {\mathrm{d} s}{m(s)}\right )\rho (t)|Y_0(t)|\,\textrm {d}t. \end{equation}
Grönwall’s inequality then implies that there is a constant
$C\gt 0$
such that
\begin{equation} |\hat {Y}_\theta (z)| \leq C|\sin \theta | \quad \text{for all $z\in [H,0]$}, \end{equation}
hence
$\lim _{\theta \to 0}\hat {Y}_\theta (z) = 0$
uniformly in
$z$
, and plugging into (A27) we see that
$\lim _{\theta \to 0}\hat {Y}_\theta '(z) = 0$
uniformly in
$z$
as well. Therefore,
\begin{equation} \lim _{\theta \to 0}\lambda _\theta = \lim _{\theta \to 0}\frac {Y_\theta '(0)}{\hat {Y}_\theta (0)} = \frac {Y_0'(0)}{Y_0(0)} = \left (\int _H^0\frac {m(0)}{m(s)}\,\textrm {d}s\right )^{-1} \gt 0, \end{equation}
whence the claim follows by arguing by continuity in (5.14).
Appendix B. Finding explicit solutions
In order to solve (5.20) with the eddy viscosity (5.23), we first determine the general solution in the upper layer
$[{z_{\star }},0]$
. To this end, we use the monotonicity of
$m$
on
$[{z_{\star }},0]$
to introduce the change of variable
\begin{equation} X(z) = \sqrt {m(z)} = \sqrt {1+\frac {M-1}{{z_{\star }}}\,z}, \quad z \in [{z_{\star }},0] \end{equation}
(Constantin & Johnson Reference Constantin and Johnson2019b ). Then, denoting
\begin{equation} W(z) = G(X), \end{equation}
the differential equation (5.20) is transformed into the Bessel equation
\begin{equation} X^2\frac {{d}^2 G}{\mathrm{d} X^2}(X) + X\frac {\mathrm{d} G}{\mathrm{d} X}(X) - \frac {8\mathrm{i} \,{z_{\star }}^2\sin \theta }{(M-1)^2} X^2 G(X) = 0, \quad X({z_{\star }})=\sqrt {M}\lt X\lt 1. \end{equation}
The general solution of (B3) is the so-called cylindrical function given by
\begin{align} G(X)=c_1J_0\left (\frac {{z_{\star }}}{M-1}\sqrt {8|\sin \theta |}\,X\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right )+c_2Y_0\left (\frac {{z_{\star }}}{M-1}\sqrt {8|\sin \theta |}\,X\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right ) \end{align}
(Polyanin & Zaitsev Reference Polyanin and Zaitsev2003), where
$c_1,c_2\in \mathbb{C}$
and
$J_0$
and
$Y_0$
are the complex-valued Bessel functions of the first and second kinds, respectively;
$J_0$
and
$Y_0$
are also known as Neumann functions (Arfken & Weber Reference Arfken and Weber2005). Therefore, given the eddy viscosity (5.23), the general solution of (5.20) (already taking into account the boundary condition (5.22)) is
\begin{align} W(z) = \begin{cases} c_1 J_0\left (\dfrac {{z_{\star }}}{M-1}\sqrt {8|\sin \theta |\left (1+\dfrac {M-1}{{z_{\star }}}\,z\right )}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right ) & \\[15pt] +\; c_2 Y_0\left (\dfrac {{z_{\star }}}{M-1}\sqrt {8|\sin \theta |\left (1+\dfrac {M-1}{{z_{\star }}}\,z\right )}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right )\!, &\quad z \in [{z_{\star }},0], \\[2em] c_3 F(z) + c_4 G(z), &\quad z\in [z_0,{z_{\star }}], \\[1em] c_5\sinh (\alpha (z-H)), &\quad z\in [H,z_0], \end{cases} \end{align}
where
$\{F,G\}$
is a basis of the solution space of the equation in the layer
$[z_0,z_\star ]$
and the coefficients
$c_1,\ldots ,c_5\in \mathbb{C}$
must be determined from the boundary conditions (5.21) and the requirement that
$W$
and
$ {\partial W}/{\partial z}$
be continuous at the interfaces
$z={z_{\star }}$
and
$z=z_0$
; moreover, we denote
\begin{equation} \alpha = (1+\textrm{i}\textrm{sign}(\theta ))\sqrt {\frac {|\sin \theta |}{\mathfrak{m}}}. \end{equation}
Enforcing the continuity of
$W$
and
$ {\partial W}/{\partial z}$
at
$z=z_0$
yields
\begin{equation} c_3 = -\frac {\zeta G(z_0) - \xi \dfrac {\partial G}{\partial z}(z_0)}{\mathscr{W}\,[F,G](z_0)}c_5 \quad \text{and} \quad c_4 = \frac {\zeta F(z_0) - \xi \dfrac {\partial F}{\partial z}(z_0)}{\mathscr{W}\,[F,G](z_0)}c_5, \end{equation}
where we abbreviated
\begin{equation} \xi = \sinh (\alpha (z_0-H)) \quad \text{and} \quad \zeta = \alpha \cosh (\alpha (z_0-H)), \end{equation}
and where
$\mathscr{W}\,[F,G]$
denotes the Wronskian of
$F$
and
$G$
\begin{equation} \mathscr{W}\,[F,G] = F\frac {\partial G}{\partial z} - \frac {\partial F}{\partial z}G, \end{equation}
which is always non-zero, since
$F$
and
$G$
are linearly independent. Without loss of generality, we may assume that
$({\partial G}/{\partial z})(z_0) \neq 0$
. Next, we require the continuity of
$W$
and
${\partial W}/{\partial z}$
at
$z=z_{\star }$
. Recalling the identities
\begin{equation} \frac {\mathrm{d} J_0(s)}{\mathrm{d} s}= -J_{1}(s) \quad \text{and} \quad \frac {\mathrm{d} Y_0(s)}{\mathrm{d} s}= -Y_{1}(s) \end{equation}
for the derivatives of the Bessel functions (Polyanin & Zaitsev Reference Polyanin and Zaitsev2003; Arfken & Weber Reference Arfken and Weber2005), a simple computation yields
\begin{equation} c_1 = \frac {\beta \hat {\delta }-\hat {\beta }\delta }{\gamma \hat {\delta } - \hat {\gamma }\delta }c_3 \quad \text{and} \quad c_2 = -\frac {\beta \hat {\gamma }-\hat {\beta }\gamma }{\gamma \hat {\delta } - \hat {\gamma }\delta }c_3, \end{equation}
with the abbreviations
\begin{equation} \begin{aligned} &\beta = F({z_{\star }}) + \kappa _{ {b}}G({z_{\star }}), \\[0.2em] &\gamma = J_0\left (\frac {{z_{\star }}}{M-1}\sqrt {8M|\sin \theta |}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right )\!, \\[0.2em] &\delta = Y_0\left (\frac {{z_{\star }}}{M-1}\sqrt {8M|\sin \theta |}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right )\!, \end{aligned} \end{equation}
and
\begin{equation} \begin{aligned} &\hat {\beta } = \frac {\partial F}{\partial z}({z_{\star }}) + \kappa _{ {b}}\frac {\partial G}{\partial z}({z_{\star }}), \\[0.2em] &\hat {\gamma } = -\sqrt {\frac {2|\sin \theta |}{M}}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}J_1\left (\frac {{z_{\star }}}{M-1}\sqrt {8M|\sin \theta |}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right )\!, \\[0.2em] &\hat {\delta } = -\sqrt {\frac {2|\sin \theta |}{M}}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}Y_1\left (\frac {{z_{\star }}}{M-1}\sqrt {8M|\sin \theta |}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right )\!, \end{aligned} \end{equation}
where
\begin{equation} \kappa _{{b}} = \frac {c_4}{c_3} = -\frac {\zeta F(z_0) - \xi \frac {\partial F}{\partial z}(z_0)}{\zeta G(z_0) - \xi \frac {\partial G}{\partial z}(z_0)}. \end{equation}
Setting
\begin{equation} \kappa = \frac {c_2}{c_1} = -\frac {\beta \hat {\gamma }-\hat {\beta }\gamma }{\beta \hat {\delta }-\hat {\beta }\delta } \quad \text{and} \quad \vartheta = \frac {{z_{\star }}}{M-1}\sqrt {8|\sin \theta |}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}, \end{equation}
with the notation of Theorem2 we find that
\begin{equation} \lambda = -\sqrt {2|\sin \theta |}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\frac {J_1(\vartheta ) + \kappa Y_1(\vartheta )}{J_0(\vartheta ) + \kappa Y_0(\vartheta )}. \end{equation}
This allows us to determine the coefficient
$c_1$
: in fact, evaluating (B5) at
$z=0$
and substituting into
\begin{equation} W(0) = \frac {\chi }{\sigma \chi + \lambda }\,W_{{w}}, \end{equation}
where
$\chi$
is the unique positive solution of (A16), we obtain
\begin{equation} c_1 = \frac {\chi }{\sigma \chi + \lambda }\boldsymbol{\cdot }\frac {W_{{w}}}{J_0(\vartheta ) + \kappa Y_0(\vartheta )}. \end{equation}
Once
$c_1$
is determined, it only remains to substitute it in the earlier identities to find all the remaining coefficients
\begin{equation} c_2 = \kappa c_1, \quad c_3 = \frac {\gamma \hat {\delta } - \hat {\gamma }\delta }{\beta \hat {\delta } - \hat {\beta }\delta }c_1, \quad c_4 = \kappa _{ {b}}c_3, \quad c_5 = -\frac {\mathscr{W}\,[F,G](z_0)}{\zeta G(z_0) - \xi \dfrac {\partial G}{\partial z}(z_0)}c_3. \end{equation}
B.1.
Constant
$m_1$
In the case
$m_1(z)=M$
, we take
$z_0 = H$
. It is then easy to see that in the layer
$[H,z_\star ]$
we have
\begin{equation} W(z) = c_3\sinh (\tilde {\alpha }(z-H)), \quad z\in [H,z_\star ], \end{equation}
where
\begin{equation} \tilde {\alpha } = (1+\textrm{i}\textrm{sign}(\theta ))\sqrt {\frac {|\sin \theta |}{M}} \quad \text{and} \quad c_3 = \frac {\gamma \hat {\delta } - \hat {\gamma }\delta }{\xi \hat {\delta } - \hat {\xi }\delta }c_1, \end{equation}
with
$\gamma ,\hat {\gamma },\delta ,\hat {\delta }$
as before and
\begin{equation} \xi = \sinh (\tilde {\alpha }(z-H)) \quad \text{and} \quad \hat {\xi } = \tilde {\alpha }\cosh (\tilde {\alpha }(z-H)). \end{equation}
The coefficients
$c_1$
and
$c_2$
are determined as discussed in the previous section, with the caveat that now
\begin{equation} \kappa = \frac {c_2}{c_1} = -\frac {\xi \hat {\gamma }-\hat {\xi }\gamma }{\xi \hat {\delta }-\hat {\xi }\delta }. \end{equation}
B.2.
Linearly decaying
$m_1$
In the case of
$m_1$
being given by (5.25), we have
\begin{equation} F(z) = J_0\left (-\dfrac {{z_{\star }}-z_0}{M-\mathfrak{m}}\sqrt {8|\sin \theta |\left (\dfrac {{z_{\star }}\mathfrak{m}-z_0 M}{{z_{\star }}-z_0} + \dfrac {M-\mathfrak{m}}{{z_{\star }}-z_0}\,z\right )}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right ) \end{equation}
and
\begin{equation} G(z) = Y_0\left (-\dfrac {{z_{\star }}-z_0}{M-\mathfrak{m}}\sqrt {8|\sin \theta |\left (\dfrac {{z_{\star }}\mathfrak{m}-z_0 M}{{z_{\star }}-z_0} + \dfrac {M-\mathfrak{m}}{{z_{\star }}-z_0}\,z\right )}\mathrm{e}^{-\textrm{i}\textrm{sign}(\theta )\frac {\pi }{4}}\right )\!. \end{equation}
B.3.
Exponentially decaying
$m_1$
With this choice (5.26) for
$m_1$
, in the layer
$[z_0,z_\star ]$
, (5.20) becomes
\begin{equation} \frac {\partial ^2 W}{\partial z^2}(z)+q\frac {\partial W}{\partial z}(z)-\frac {2\mathrm{i}\sin \theta }{M}\mathrm{e}^{-q(z-{z_{\star }})}W(z)=0, \quad z_0 \lt z \lt {z_{\star }}. \end{equation}
Let us define the function
$H$
through
\begin{equation} W(z) = H(z)\mathrm{e}^{-\frac {q}{2}(z-{z_{\star }})}. \end{equation}
Then
\begin{equation} \frac {\partial W}{\partial z}(z) = \mathrm{e}^{-\frac {q}{2}(z-{z_{\star }})}\left (H'(z)-\frac {q}{2}H(z)\right ) \end{equation}
and
\begin{equation} \frac {\partial ^2 W}{\partial z^2}(z) = \mathrm{e}^{-\frac {q}{2}(z-{z_{\star }})}\left (\frac {\partial ^2 H}{\partial z^2}(z) - q\frac {\partial H}{\partial z}(z) +\frac {q^2}{4}H(z)\right )\!, \end{equation}
and (B26) transforms into
\begin{equation} \frac {\partial ^2 H}{\partial z^2}(z) - \left (\frac {2\mathrm{i}\sin \theta }{M}\mathrm{e}^{-q(z-{z_{\star }})} + \frac {q^2}{4}\right )H(z) = 0, \quad z_0 \lt z \lt 0. \end{equation}
Introducing the change of variables
\begin{equation} X = \frac {2}{q}\sqrt {\frac {2\mathrm{i} \sin \theta }{M}} \mathrm{e}^{-\frac {q}{2}(z-{z_{\star }})}, \end{equation}
for which we have
\begin{equation} \frac {{d}}{\mathrm{d} z} = -\frac {q}{2}X\frac {{d}}{\mathrm{d} X} \quad \text{and} \quad \frac {{d}^2}{\mathrm{d} z^2} = \frac {q^2}{4}\left (X^2\frac {{d}^2}{\mathrm{d} X^2} + X\frac {{d}}{\mathrm{d} X}\right )\!, \end{equation}
and denoting
\begin{equation} H(z) = G(X), \end{equation}
then (B30) transforms into the equation
\begin{equation} X^2\frac {{d}^2 G}{\mathrm{d} X^2}(X) + X\frac {\mathrm{d} G}{\mathrm{d} X}(X) - (X^2+1)G(X) = 0, \quad \frac {2}{q}\sqrt {\frac {2\mathrm{i} \sin \theta }{M}} \lt X \lt \frac {2}{q}\sqrt {\frac {2\mathrm{i} \sin \theta }{\mathfrak{m}}}, \end{equation}
whose general solution is
\begin{equation} G(X) = c_3 I_1(X) + c_4 K_1(X), \quad \frac {2}{q}\sqrt {\frac {2\mathrm{i} \sin \theta }{M}} \lt X \lt \frac {2}{q}\sqrt {\frac {2\mathrm{i} \sin \theta }{\mathfrak{m}}}, \end{equation}
where
$c_3,c_4\in \mathbb{C}$
and
$I_1$
and
$K_1$
are the modified Bessel functions of the first and second kinds, respectively, both of order 1. Hence,
\begin{equation} F(z) = \mathrm{e}^{-\frac {q}{2}(z-{z_{\star }})}I_1\left (\frac {2}{q}\sqrt {\frac {2\mathrm{i} \sin \theta }{M}} \mathrm{e}^{-\frac {q}{2}(z-{z_{\star }})}\right ) \quad \text{and} \quad G(z) = \mathrm{e}^{-\frac {q}{2}(z-{z_{\star }})}K_1\left (\frac {2}{q}\sqrt {\frac {2\mathrm{i} \sin \theta }{M}} \mathrm{e}^{-\frac {q}{2}(z-{z_{\star }})}\right )\!. \end{equation}
Recalling that
\begin{equation} \frac {\mathrm{d} I_1(\mathrm{t})}{\mathrm{dt}} = \frac {1}{2}( I_0(\mathrm{t})+ I_2(\mathrm{t})) \quad \text{and} \quad \frac {\mathrm{d} K_1(\mathrm{t})}{\mathrm{dt}} = -\frac {1}{2}( K_0(\mathrm{t})+ K_2(\mathrm{t})), \end{equation}
a straightforward computation as described above yields the coefficients
$c_1,\ldots ,c_5$
.
Appendix C. The problem at first order in the expansion with respect to
$\boldsymbol{\mathscr{R}}$
In this appendix, we derive the equations for the first-order correction of the leading-order dynamics. Once more setting the coefficients to zero in the asymptotic expansions (3.11), (3.12), (3.13), (3.14) and (3.15), we obtain the problem
\begin{equation} \left .\begin{aligned} &\frac {1}{\rho }\frac {\partial }{\partial z}\left (m\frac {\partial \hat {u}}{\partial z}\right ) + 2\hat {v}\sin \theta \\ &\quad \quad = \frac {1}{\rho \cos \theta }\frac {\partial \hat {p}}{\partial \varphi } -\left (\frac {\tilde {u}}{\cos \theta }\frac {\partial \tilde {u}}{\partial \varphi } + \tilde {v}\frac {\partial \tilde {u}}{\partial \theta } + \tilde {w}\frac {\partial \tilde {u}}{\partial z}\right ) + \tilde {u}\tilde {v}\tan \theta \\ &\frac {1}{\rho }\frac {\partial }{\partial z}\left (m\frac {\partial \hat {v}}{\partial z} \right ) -2\hat {u}\sin \theta \\ &\quad \quad = \frac {1}{\rho }\frac {\partial \hat {p}}{\partial \theta } - \left (\frac {\tilde {u}}{\cos \theta }\frac {\partial \tilde {v}}{\partial \varphi } + \tilde {v}\frac {\partial \tilde {v}}{\partial \theta } + \tilde {w}\frac {\partial \tilde {v}}{\partial z}\right ) - \tilde {u}^2\tan \theta \\ &\frac {\partial \hat {p}}{\partial z} = 0 \\ &\frac {\partial \hat {u}}{\partial \varphi } + \frac {\partial }{\partial \theta } \left (\hat {v} \cos \theta \right ) + \cos \theta \frac {\partial \hat {w}}{\partial z} = 0 \end{aligned}\;\right \} \quad \text{in} \quad \{H \lt z \lt 0\}, \end{equation}
with dynamic boundary condition
\begin{equation} \hat {p} = 0 \quad \text{on} \quad \{z=0\} \end{equation}
and upper boundary condition due to the wind stress given by
\begin{equation} \begin{pmatrix}\dfrac {\partial \hat {u}}{\partial z}\\[1em] \dfrac {\partial \hat {v}}{\partial z}\end{pmatrix} = -\sigma \left (\sqrt {(u_{{a}}-\sigma \tilde {u})^2+(v_{{a}}-\sigma \tilde {v})^2}\begin{pmatrix} \hat {u} \\[0.2em] \hat {v} \end{pmatrix} + \frac {\hat {u}(u_{{a}} - \sigma \tilde {u}) + \hat {v}(v_{{a}} - \sigma \tilde {v})}{\sqrt {(u_{{a}}-\sigma \tilde {u})^2+(v_{\textrm { a}}-\sigma \tilde {v})^2}}\begin{pmatrix}u_{{a}}-\sigma \tilde {u} \\[0.2em] v_{{a}}-\sigma \tilde {v}\end{pmatrix}\right ) \quad \text{on} \quad \{z = 0\}. \end{equation}
The equation
$ {\partial \hat {p}}/{\partial z} = 0$
, together with (C2), implies that
$\hat {p}=0$
throughout. Therefore, we are left with the problem
\begin{equation} \left .\begin{aligned} \frac {\partial }{\partial z}\left (m\frac {\partial \hat {u}}{\partial z}\right ) + 2\rho \hat {v}\sin \theta &= \mathscr{F}(\tilde {u},\tilde {v}) \\[0.2em] \frac {\partial }{\partial z}\left (m\frac {\partial \hat {v}}{\partial z}\right ) + 2\rho \hat {u}\sin \theta &= \mathscr{G}(\tilde {u},\tilde {v}) \end{aligned}\right \} \quad \text{in} \quad \{H \lt z \lt 0\}, \end{equation}
where
\begin{equation} \mathscr{F}(\tilde {u},\tilde {v}) = - \rho \left (\frac {\tilde {u}}{\cos \theta }\frac {\partial \tilde {u}}{\partial \varphi } + \tilde {v}\frac {\partial \tilde {u}}{\partial \theta } + \tilde {w}\frac {\partial \tilde {u}}{\partial z}\right ) + \rho \tilde {u}\tilde {v}\tan \theta \end{equation}
and
\begin{equation} \mathscr{G}(\tilde {u},\tilde {v}) = - \rho \left (\frac {\tilde {u}}{\cos \theta }\frac {\partial \tilde {v}}{\partial \varphi } + \tilde {v}\frac {\partial \tilde {v}}{\partial \theta } + \tilde {w}\frac {\partial \tilde {v}}{\partial z}\right ) - \rho \tilde {u}^2\tan \theta . \end{equation}
The leading-order velocity
$(\tilde {u},\tilde {v})$
enters the equations only through the inhomogeneity, given by
$\mathscr{F}$
and
$\mathscr{G}$
, and the upper boundary condition. Once the solution for
$(\hat {u},\hat {v})$
has been obtained from these equations, from the continuity equation we can recover
$\hat {w}$
via
\begin{equation} \frac {\partial \hat {w}}{\partial z} = -\frac {1}{\cos \theta }\frac {\partial \hat {u}}{\partial \varphi } - \frac {\partial \hat {v}}{\partial \theta } +\tan \theta \, \hat {v}. \end{equation}
As with the leading-order problem, it is convenient to introduce the complex notation
\begin{equation} \hat {W} = \hat {u} + \mathrm{i}\hat {v} \quad \text{and} \quad \mathscr{I} = \mathscr{F} + \mathrm{i}\mathscr{G}; \end{equation}
then, (C4) becomes
\begin{equation} \frac {\partial }{\partial z}\left (m\frac {\partial \hat {W}}{\partial z}\right ) - 2\mathrm{i}\rho \hat {W}\sin \theta = \mathscr{I}. \end{equation}
Moreover, with this notation, and using the identity
$\mathrm{Re}(\alpha ) = ({1}/{2})(\alpha + \bar {\alpha })$
for
$\alpha \in \mathbb{C}$
, we can rewrite the boundary condition (C3) as
\begin{equation} \begin{aligned} \frac {\partial \hat {W}}{\partial z} &= -\sigma \left (|W_{{w}} - \sigma W|\hat {W} + \frac {W_{{w}} - \sigma W}{|W_{{w}} - \sigma W|}\mathrm{Re}\left \{(W_{{w}} - \sigma W\overline {\hat {W}}\right \}\right ) & & \\ &= -\frac {\mu }{2}\left (3|W_{{w}} - \sigma W|\hat {W} + \frac {(W_{{w}} - \sigma W^2}{|W_{{w}} - \sigma W|}\overline {\hat {W}}\right ) & & \text{on $\{z=0\}$}, \end{aligned} \end{equation}
where
$W_{{w}}$
and
$W$
are as in (5.6), and the bar denotes complex conjugation. Thus, also taking (3.7) into account, we are left with the problem
\begin{equation} \left \{\begin{aligned} & \frac {\partial }{\partial z}\left (m\frac {\partial \hat {W}}{\partial z}\right ) - 2\mathrm{i}\rho \hat {W}\sin \theta = \mathscr{I} & &\quad \text{in} \quad \{H\lt z\lt 0\}, \\[0.2em] & \frac {\partial \hat {W}}{\partial z} = -\frac {\mu }{2}\left (3|W_{{w}} - \sigma W|\hat {W} + \frac {(W_{{w}} - \sigma W^2}{|W_{{w}} - \sigma W|}\overline {\hat {W}}\right ) & &\quad \text{on} \quad \{z=0\}, \\[0.2em] & \hat {W} = 0 & &\quad \text{on} \quad \{z=H\}. \end{aligned}\right . \end{equation}
In order to solve (C11), the standard solution method is the variation of constants. First, we look for a basis
$\{\hat {W}_1,\hat {W}_2\}$
of the solution space of the homogeneous equation
\begin{equation} \frac {\partial }{\partial z}\left (m\frac {\partial \hat {W}}{\partial z}\right ) - 2\mathrm{i}\rho \hat {W}\sin \theta = 0 \quad \text{in} \quad \{H\lt z\lt 0\}. \end{equation}
Then, a particular solution
$\hat {W}_{ {p}}$
of the inhomogeneous equation (C9) is given by
\begin{equation} \hat {W}_{{p}} = w_1\hat {W}_1 + w_2\hat {W}_2; \end{equation}
here, the coefficients
$w_1$
and
$w_2$
are functions that are determined by integrating
\begin{equation} \frac {\partial w_1}{\partial z} = \frac {\hat {W}_2\mathscr{I}}{\mathscr{W}\,[\hat {W}_1,\hat {W}_2]} \quad \text{and} \quad \frac {\partial w_2}{\partial z} = -\frac {\hat {W}_1\mathscr{I}}{\mathscr{W}\,[\hat {W}_1,\hat {W}_2]}, \end{equation}
where, as before,
\begin{equation} \mathscr{W}\,[\hat {W}_1,\hat {W}_2] = \hat {W}_1\frac {\partial \hat {W}_2}{\partial z} -\frac {\partial \hat {W}_1}{\partial z}\hat {W}_2 \end{equation}
denotes the Wronskian of
$\hat {W}_1$
and
$\hat {W}_2$
(which is always non-zero, because
$\hat {W}_1$
and
$\hat {W}_2$
are linearly independent). Without loss of generality, we can choose
$w_1(H) = w_2(H) = 0$
. The general solution of (C9) is then given by
\begin{equation} \hat {W} = C_1\hat {W}_1 + C_2\hat {W}_2 + \hat {W}_{\textrm {p}}, \end{equation}
where
$C_1,C_2\in \mathbb{C}$
are arbitrary constants, to be determined through the boundary conditions. In fact, assuming without loss of generality that
$\hat {W}_2(H) \neq 0$
, the third equation in (C11), together with
$\hat {W}_{{p}}(H)=0$
, yields
\begin{equation} C_2 = -\frac {\hat {W}_1(H)}{\hat {W}_2(H)}C_1, \end{equation}
hence
\begin{equation} \hat {W}(z) = C_1\left (\hat {W}_1(z) - \frac {\hat {W}_1(H)}{\hat {W}_2(H)}\hat {W}_2(z)\right ) + \hat {W}_{\textrm {p}}(z). \end{equation}
Then, plugging into the second equation of (C11), we obtain the following equation for
$C_1$
:
\begin{equation} \alpha C_1 + \beta \overline {C_1} = \gamma , \end{equation}
where, as before, the bar denotes complex conjugation, and
\begin{equation} \begin{aligned} \alpha &= \hat {W}_1'(0) - \frac {\hat {W}_1(H)}{\hat {W}_2(H)}\hat {W}_2'(0) + \frac {3}{2}\sigma |W_{{w}} - \sigma W(0)|\left (\hat {W}_1(0) - \frac {\hat {W}_1(H)}{\hat {W}_2(H)}\hat {W}_2(0)\right )\!, \\ \beta &= \frac {1}{2}\sigma \frac {(W_{{w}} - \sigma W(0))^2}{|W_{{w}} - \sigma W(0)|}\left (\overline {\hat {W}_1(0)} - \frac {\overline {\hat {W}_1(H)}}{\overline {\hat {W}_2(H)}}\overline {\hat {W}_2(0)}\right )\!, \\ \gamma &= -\frac {1}{2}\sigma \biggl [3|W_{{w}} - \sigma W(0)|\bigl (w_1(0)\hat {W}_1(0) + w_2(0)\hat {W}_2(0)\bigr )\biggr . \\ &\hspace {1.2cm} \biggl .+ \frac {(W_{{w}} - \sigma W(0))^2}{|W_{{w}} - \sigma W(0)|}\left (\overline {w_1(0)}\overline {\hat {W}_1(0)} + \overline {w_2(0)}\overline {\hat {W}_2(0)}\right )\biggr ] \\ &\quad - w_1(0)\hat {W}_1'(0) - w_2(0)\hat {W}_2'(0). \end{aligned} \end{equation}
Equation (C19) has the unique solution
\begin{equation} C_1 = \frac {\bar {\alpha }\gamma - \beta \bar {\gamma }}{|\alpha | - |\beta |}, \end{equation}
as long as
$|\alpha | \neq |\beta |$
. This condition can be reformulated as follows. Denote
\begin{equation} Y(z) = \hat {W}_1(z) - \frac {\hat {W}_1(H)}{\hat {W}_2(H)}\hat {W}_2(z); \end{equation}
then, assuming without loss of generality that
$Y'(H)=1$
, the function
$Y$
is the solution of (A1). Thus, writing
\begin{equation} \lambda = \frac {Y'(0)}{Y(0)} \end{equation}
and recalling that
$\chi = |W_{{w}} - \sigma W(0)|$
is given by the unique positive solution of (A16), as in the proof of Theorem1, we see that
$|\alpha | \neq |\beta |$
is equivalent to
\begin{equation} \left |\lambda + \frac {3}{2}\sigma \chi \right | \neq \frac {1}{2}\sigma \chi . \end{equation}
Since
\begin{equation} \left |\lambda + \frac {3}{2}\sigma \chi \right | \ge |\lambda | - \frac {3}{2}\sigma \chi \end{equation}
and, as we argued in the proof of Theorem2,
\begin{equation} \sigma \chi \le \frac {-\mathrm{Re}(\lambda ) + \sqrt {\mathrm{Re}(\lambda )^2 + 4\sigma |\lambda W_{{w}}|}}{2} \le \sigma |W_{{w}}|\frac {|\lambda |}{\mathrm{Re}(\lambda )}, \end{equation}
where we used the inequality
$\sqrt {1+x} \le 1+({1}/{2})x$
for
$x \ge -1$
, combining (C25) and (C26) we have
\begin{equation} \left |\lambda + \frac {3}{2}\sigma \chi \right | \ge \left (\frac {\mathrm{Re}(\lambda )}{|W_{{w}}|} - \frac {3}{2}\sigma \right )\chi . \end{equation}
Therefore, the validity of (C24) is guaranteed if
\begin{equation} \mathrm{Re}(\lambda ) \gt 2\sigma |W_{{w}}|. \end{equation}
Of course, this in fact boils down to a condition on the eddy viscosity and the wind’s intensity, since
$\lambda$
depends only on the former. However, finding explicit bounds on
$\textrm { Re}(\lambda )$
in terms of a general
$m$
is difficult, so in practice, this condition would have to be checked case by case. In this place, let us only note that the right-hand side of (C28) is very small, of order of magnitude
$10^{-2}$
, and in the explicit cases that we considered earlier, the condition is satisfied by a large margin.
It is also possible to provide estimates for the first-order correction of the horizontal velocity field given by (C11). To this end, we recall the notion of the logarithmic matrix norm.
Definition 1.
For
$n\in \mathbb{N}$
, let
$\|\boldsymbol{\cdot }\|_*$
denote a matrix norm on
$\mathbb{C}^{n\times n}$
induced by a vector norm
$\|\boldsymbol{\cdot }\|$
on
$\mathbb{C}^n$
. For a matrix
$A\in \mathbb{C}^{n\times n}$
, we define its logarithmic matrix norm
$\mu _*[A]$
associated with
$\|\boldsymbol{\cdot }\|_*$
as
\begin{equation} \mu _*[A] = \lim _{h\searrow 0}\frac {\|\mathbb{I}+hA\|_*-1}{h}, \end{equation}
where
$\mathbb{I}$
denotes the identity matrix in
$\mathbb{C}^{n\times n}$
.
It is easy to show that
$\mu _*[A]$
is well defined for all
$A\in \mathbb{C}^{n\times n}$
. For more properties of the logarithmic matrix norm, we refer the reader to Söderlind (Reference Söderlind2024). These logarithmic matrix norms were first applied to the classical Ekman problem by Marynets (Reference Marynets2022). In fact, they can be especially useful, among other things, in determining bounds on solutions to ODEs. The following result, which will be useful to us, can be found in § 6.6 of Söderlind (Reference Söderlind2024).
Lemma 1.
For
$T\gt 0$
, let
$A:[0,T)\to \mathbb{C}^{n\times n}$
and
$f:[0,T)\to \mathbb{C}^n$
both be continuous. Suppose that
$u\in C^1((0,T),\mathbb{C}^n)\cap C([0,T),\mathbb{C}^n)$
solves the inhomogeneous linear initial value problem
\begin{equation} u'(t) = A(t)u(t) + f(t), \quad u(0)=u_0. \end{equation}
Denote
\begin{equation} L(t) = \int _0^t\mu _*[A(\tau )]\,\textrm {d}\tau . \end{equation}
Then
\begin{equation} \|u(t)\| \leq \mathrm{e}^{L(t)}\|u_0\| + \mathrm{e}^{L(t)}\int _0^t\mathrm{e}^{-L(s)}\|f(s)\|\,\textrm {d}s. \end{equation}
In order to apply Lemma1, we rewrite the first equation of (C11) as
\begin{equation} \frac {\partial x}{\partial z}(z) = A(z)x(z) + f(z), \quad z\in (H,0), \end{equation}
where
\begin{align} x = \begin{pmatrix} \hat {W} \\[0.2em] \dfrac {\partial \hat {W}}{\partial z} \end{pmatrix}\!, \quad A = \begin{pmatrix} 0 & 1 \\[0.2em] \dfrac {2\mathrm{i}\sin \theta \,\rho }{m} & -\dfrac {m'}{m} \end{pmatrix}\!, \quad \text{and} \quad f = \begin{pmatrix} 0 \\[0.2em] \dfrac {\mathscr{I}}{m} \end{pmatrix}\!. \end{align}
Choosing the norm
$\|x\|_1 = |x_1| + |x_2|$
and denoting
$A = (a_{jk})_{1\le j,k\le 2}$
, the corresponding logarithmic matrix norms can be computed explicitly as
\begin{equation} \mu _1[A(z)] = \max _k\biggl \{\mathrm{Re}(a_{kk}(z)) + \sum _{j\neq k}|a_{jk}(z)|\biggr \} = \max \left \{2|\sin \theta |,\,1 - \frac {m'(z)}{m(z)}\right \} \end{equation}
(Söderlind Reference Söderlind2024). By Lemma1 and (C35), we have
\begin{align} \bigl |\hat {W}(z)\bigr | + \left |\frac {\partial \hat {W}}{\partial z}(z)\right | \leq \left |\frac {\partial \hat {W}}{\partial z}(H)\right |\exp \left (\int _{H}^z\mu _1[A(s)]\,\textrm {d}s\right ) + \int _{H}^z\exp \left (\int _t^z\mu _1[A(s)]\,\textrm {d}s\right )\frac {|\mathscr{I}(t)|}{m(t)}\,\textrm {d}t. \end{align}
Note that this bound depends only on the eddy viscosity
$m$
and, through
$\mathscr{I}$
and
$( {\partial \hat {W}}/{\partial z})(H)$
, the leading-order solution
$W$
found in § 5.1.
Appendix D. Ekman spirals
In this appendix, we include some three-dimensional plots of the Ekman spirals for the three profiles of the eddy viscosity considered in § 5.2. More precisely, figure 11 refers to the case of constant
$m_1$
, whereas figure 12 shows the case of a linearly decaying
$m_1$
and, lastly, figure 13 displays the case of an exponentially decaying
$m_1$
.
The case of constant
$m_1$
for
$|W_{{w}}|=1$
,
$\sigma = 0.014$
(corresponding to
$U_{{a}}' = 10\,\mathrm{m\,s^{-1}}$
) and
$H=-80$
: (a)
$M = 10$
,
${z_{\star }} = -1$
, (b)
$M = 50$
,
${z_{\star }} = -1$
.
Figure 11. Long description
Panel A: A 3D line graph with axes labeled as z, ue, and ve. The graph shows multiple lines in different colors (red, blue, black) representing different data series. The values on the axes range from -50 to 0 for z, -0.20 to 0.25 for ue, and -0.10 to 0.15 for ve. The lines form a 3D shape indicating the relationship between these variables for M = 10 and z* = -1. Panel B: Another 3D line graph with similar axes labels and color schemes. The values on the axes range from -50 to 0 for z, -0.08 to 0.14 for ue, and -0.04 to 0.10 for ve. This graph represents the relationship for M = 50 and z* = -1. Both graphs illustrate how the variables interact in different scenarios.
The case of linearly decaying
$m_1$
for
$|W_{{w}}|=1$
,
$\sigma = 0.014$
,
$H=-80$
and
$\mathfrak{m}=0.01$
: (a)
$M=10$
,
$z_0=-5$
,
$z_{\star }=-1$
, (b)
$M=10$
,
$z_0=-20$
,
$z_{\star }=-4$
, (c)
$M=50$
,
$z_0=-5$
,
$z_{\star }=-1$
, (d)
$M=50$
,
$z_0=-20$
,
$z_{\star }=-4$
.
The case of exponentially decaying
$m_1$
for
$|W_{{w}}|=1$
,
$\sigma = 0.014$
,
$H=-80$
and
$\mathfrak{m}=0.01$
: (a)
$M=10$
,
$z_0=-5$
,
$z_{\star }=-1$
, (b)
$M=10$
,
$z_0=-20$
,
$z_{\star }=-4$
, (c)
$M=50$
,
$z_0=-5$
,
$z_{\star }=-1$
, (d)
$M=50$
,
$z_0=-20$
,
$z_{\star }=-4$
.
Figure 13. Long description
Panel A: A 3D line graph with axes labeled as z, ue, and ve. The graph shows three curves in different colors: red, blue, and black. The parameters are M = 10, z0 = -5, and z* = -1. The curves represent different relationships between the variables. Panel B: A 3D line graph with axes labeled as z, ue, and ve. The graph shows three curves in different colors: red, blue, and black. The parameters are M = 10, z0 = -20, and z* = -4. The curves represent different relationships between the variables. Panel C: A 3D line graph with axes labeled as z, ue, and ve. The graph shows three curves in different colors: red, blue, and black. The parameters are M = 50, z0 = -5, and z* = -1. The curves represent different relationships between the variables. Panel D: A 3D line graph with axes labeled as z, ue, and ve. The graph shows three curves in different colors: red, blue, and black. The parameters are M = 50, z0 = -20, and z* = -4. The curves represent different relationships between the variables.
ms−1
W=u~e+iv~e
m1
|Ww|=1
σ=0.014
Ua′=10ms−1
θ=45∘
H=−80
z=−30
M=10
z⋆=−1
M=50
z⋆=−1
m1
|Ww|=1
σ=0.014
H=−80
M=10
M=50
W=u~e+iv~e
m1
|Ww|=1
σ=0.014
H=−80
m=0.01
z=−21
M=10
z0=−5
z⋆=−1
M=10
z0=−20
z⋆=−4
M=50
z0=−5
z⋆=−1
M=50
z0=−20
z⋆=−4
m1
|Ww|=1
σ=0.014
H=−80
m=0.01
M=10
z0=−5
M=10
z0=−20
M=50
z0=−5
M=50
z0=−20
W=u~e+iv~e
m1
|Ww|=1
σ=0.014
z0=−5
H=−80
m=0.01
z=−12
M=10
z0=−5
z⋆=−1
M=10
z0=−20
z⋆=−4
M=50
z0=−5
z⋆=−1
M=50
z0=−20
z⋆=−4
m1
|Ww|=1
σ=0.014
H=−80
m=0.01
M=10
z0=−5
M=10
z0=−20
M=50
z0=−5
M=50
z0=−20
m1
|Ww|=1
σ=0.014
z0=−5
H=−80
m=0.01
z=−5
M=10
z⋆=−0.1
M=10
z⋆=−1
M=50
z⋆=−0.1
M=50
z⋆=−1
m1
|Ww|=1
σ=0.014
Ua′=10ms−1
H=−80
M=10
z⋆=−1
M=50
z⋆=−1
m1
|Ww|=1
σ=0.014
H=−80
m=0.01
M=10
z0=−5
z⋆=−1
M=10
z0=−20
z⋆=−4
M=50
z0=−5
z⋆=−1
M=50
z0=−20
z⋆=−4
m1
|Ww|=1
σ=0.014
H=−80
m=0.01
M=10
z0=−5
z⋆=−1
M=10
z0=−20
z⋆=−4
M=50
z0=−5
z⋆=−1
M=50
z0=−20
z⋆=−4
Open access
