1 Introduction

Wind-driven currents are one of the main aspects of oceanic circulation, which is largely responsible for the heat transport between different geographical regions and has a crucial impact on the Earth’s climate [15]. At non-equatorial latitudes, the flow below a depth of about 100 metres is approximately in geostrophic balance, that is, the force balance occurs between the Coriolis force and horizontal pressure gradients. The frictional effects of surface winds are confined to a relatively thin boundary layer, called the Ekman layer, where the momentum of the wind is carried down to the interior of the water flow by the wind-induced turbulence of the water. Thus, in this surface layer we find currents (superimposed to the underlying geostrophic motion) where the force balance is between the Coriolis force and the frictional forces due to the wind (see [15, 18]). It is observed in practice that such wind-induced currents do not follow the direction of the generating wind, but instead are deflected to the right of the wind direction in the Northern Hemisphere (to the left in the Southern Hemisphere), the amplitude of the observed deflection angle being generally within the range 10\(^\circ \)–75\(^\circ \) (cf. the data in [17] or [20]).

On the other hand, the ocean dynamics near the Equator differs significantly from that at higher latitudes. In fact, since the Coriolis force vanishes along the Equator, geostrophic effects are negligible there and nonlinear effects arise at leading order; moreover, the Coriolis force changes sign across the Equator and so acts as a waveguide that gives rise to azimuthal (i.e. longitudinal) flow propagation (see e.g. [3, 4, 6, 12]). Therefore, substantially different approaches are required for the description of non-equatorial and equatorial flows. This paper is devoted to the former aspect of oceanic wind-driven currents; for the latter case, we refer the reader to the references above.

The Ekman layer is named after the Swedish scientist who, in 1905, was the first to write down and investigate a model that should explain the observations of the surface deflection angle that had been made at that time. In his classical work [10], the eddy viscosity (a quantity introduced as a measure of the frictional effects in the boundary layer—see the next section) is set as a constant, yielding a solution that, at the surface, is directed at an angle of 45\(^\circ \) to the right (left) of the wind direction in the Northern (Southern) Hemisphere. The source of the inaccuracy of the predicted surface deflection angle has been identified in the simplified assumption of constant eddy viscosity, and over the last few decades many attempts have been made at finding solutions for more general depth-dependent eddy viscosity profiles; many of the known explicit solutions for particular cases can be found for instance in [5, 7, 8, 11, 14, 16]. Nevertheless, finding explicit solutions can be a very hard task even for rather simple eddy viscosity profiles, so alternative ways to tackle the problem have been sought. In [19], an application of the WKB approach yields an asymptotic approximation of the solution to the general Ekman model; in [1, 2], a formula for the deviation of the deflection angle from the reference value of 45\(^\circ \) for a given small perturbation of a constant eddy viscosity is derived; in [9], the problem of a two-layer ocean with constant eddy viscosity in the lower layer and a general continuous eddy viscosity in the upper layer is transferred to solving a Riccati-type equation in the upper layer, from whose solution the surface deflection angle for the original problem can be determined straightforwardly. In the present article we combine the latter two methods: using the reformulation of the problem as a Riccati-type equation, we perform a perturbation analysis to derive a formula for the change of the surface deflection angle caused by a perturbation of a general continuous eddy viscosity profile in the upper layer.

2 The Governing Equations

The steady-state equations for non-equatorial Ekman-type flows in the f-plane approximation for geophysical flows (cf. [15, 18]) read as follows (see [2, 8]):

$$\begin{aligned} -fv&= \frac{\mathrm {d}}{\mathrm {d}Z}\left( \nu \frac{\mathrm {d}u}{\mathrm {d}Z}\right) \\ fu&= \frac{\mathrm {d}}{\mathrm {d}Z}\left( \nu \frac{\mathrm {d}v}{\mathrm {d}Z}\right) , \end{aligned}$$

where (uv) are the horizontal Cartesian components of the flow in the f-plane, \(Z\in (-\infty ,0)\) is the depth below the mean surface level (located at \(Z=0\)), \(\nu \) is the (depth-dependent) eddy viscosity, and \(f = 2\Omega \sin (\theta )\) is the Coriolis parameter at the fixed latitude \(\theta \in (-\pi /2,0)\cup (0,\pi /2)\), with the angular speed of rotation of the Earth \(\Omega = 7.29\times 10^5\,\)rad/s; in the Northern Hemisphere we have \(f>0\), whereas \(f<0\) in the Southern Hemisphere. In the following we choose to work in the Northern Hemisphere; the results can be readily adapted to the other case. Introducing the complex notation \(U = u+\mathrm {i}v\), the equations above can be rewritten as

$$\begin{aligned} \mathrm {i}f U = \frac{\mathrm {d}}{\mathrm {d}Z}\left( \nu \frac{\mathrm {d}U}{\mathrm {d}Z}\right) \quad \text {in } (-\infty ,0). \end{aligned}$$
(2.1)

The surface boundary condition is

$$\begin{aligned} \rho \nu (0)\frac{\mathrm {d}U}{\mathrm {d}Z}(0) = \tau _1+\mathrm {i}\tau _2, \end{aligned}$$
(2.2)

where \(\rho \) is the (constant) density of the water and \(\tau _1+\mathrm {i}\tau _2\) is the surface wind stress, which we assume to be constant as well. In fact, given the fact that in the f-plane approximation we are dealing with a planar geometry and constant Coriolis parameter f, we can rotate the horizontal coordinate system so that the real axis becomes aligned with the wind direction: this means that, without loss of generality, we may assume \(\tau _2=0\) and \(\tau _1>0\). Moreover, we require

$$\begin{aligned} U(Z) \rightarrow 0 \quad \text {as } Z\rightarrow -\infty . \end{aligned}$$
(2.3)

The problem (2.1)–(2.3) can be non-dimensionalized by setting

$$\begin{aligned} K = \frac{f\rho }{\tau _1}\nu , \quad \psi = \frac{K(0)}{\sqrt{2\tau _1/\rho }}U, \quad z = \frac{f}{\sqrt{2\tau _1/\rho }}Z.\end{aligned}$$

From this we obtain the non-dimensional problem

$$\begin{aligned} (K(z)\psi '(z))' -2\mathrm {i}\psi (z) = 0&\quad \text {for } z<0 \end{aligned}$$
(2.4a)
$$\begin{aligned} \psi ' = 1&\quad \text {on } z=0 \end{aligned}$$
(2.4b)
$$\begin{aligned} \psi (z)\rightarrow 0&\quad \text {as } z\rightarrow -\infty , \end{aligned}$$
(2.4c)

where the prime denotes a derivative with respect to z. Let us remark that the surface deflection angle, which is the main objective of this paper, is the argument of the complex vector \(\psi (0)\).

Now, let us suppose that, for the continuous \(K_0:(-\infty ,0]\rightarrow (0,\infty )\) with \(K_0(z)=1\) for \(z\le -h\) (for some non-dimensional depth \(-h<0\)), the solution \(\psi _0\) to

$$\begin{aligned} (K_0(z)\psi _0'(z))' -2\mathrm {i}\psi _0(z) = 0&\quad \text {for } z<0 \\ \psi _0' = 1&\quad \text {on } z=0 \\ \psi _0(z)\rightarrow 0&\quad \text {as } z\rightarrow -\infty \end{aligned}$$

is known. Henceforth we consider a K of the form

$$\begin{aligned} K(z) = K_0(z) + \varepsilon K_1(z) + o(\varepsilon ) >0, \quad z \le 0, \end{aligned}$$
(2.5)

for small \(\varepsilon >0\), uniformly in z: explicitly, this means that for each \(C>0\) there exists \(\varepsilon _0\) such that for all \(\varepsilon \le \varepsilon _0\) and all \(z\le 0\) we have

$$\begin{aligned} \vert K(z) - [K_0(z) + \varepsilon K_1(z)]\vert \le C\varepsilon .\end{aligned}$$

Moreover, we assume that \(K_1\) is continuous on \((-\infty ,0]\) and

$$\begin{aligned} K_1(z)=0, \quad z\le -h.\end{aligned}$$

Notice that the requirement of a constant K below some depth is physically reasonable (cf. [9]): in fact, the turbulence is confined to the thin Ekman boundary layer, thus the eddy viscosity at greater depth coincides with that of sea water, which can be normalized to obtain

$$\begin{aligned}K(z)=1, \quad z\le -h.\end{aligned}$$

3 Perturbation Analysis

In [9] it is shown that, if q solves the Riccati equation

$$\begin{aligned} q'(z)+\frac{q(z)^2}{K(z)}&= 2\mathrm {i}, \quad z\in (-h,0) \end{aligned}$$
(3.1a)
$$\begin{aligned} q(-h)&= 1+\mathrm {i}, \end{aligned}$$
(3.1b)

then for the restriction of the solution \(\psi \) of (2.4) on \((-h,0)\) we have that \(\psi (0)=K(0)/q(0)\),

$$\begin{aligned} \arg (\psi (0)) = -\arg (q(0)) \end{aligned}$$
(3.2)

and

$$\begin{aligned} \psi (z) = \frac{K(0)}{q(0)}\exp \left( -\int \limits _z^0\frac{q(y)}{K(y)}\,\mathrm {d}y\right) , \quad z\in (-h,0). \end{aligned}$$
(3.3)

Following (2.5), let us make the ansatz

$$\begin{aligned} q(z) = q_0(z) + \varepsilon q_1(z) + o(\varepsilon ) \end{aligned}$$
(3.4)

(uniformly in z) for the solution to (3.1), where \(q_0\) solves

$$\begin{aligned} q_0'(z) + \frac{q_0(z)^2}{K_0(z)}&= 2\mathrm {i}, \quad z\in (-h,0) \end{aligned}$$
(3.5a)
$$\begin{aligned} q_0(-h)&= 1+\mathrm {i}, \end{aligned}$$
(3.5b)

and \(q_1\) is to be determined. Inserting (2.5) and (3.4) into (3.1a), using (3.5a), dividing by \(\varepsilon \), and letting \(\varepsilon \rightarrow 0\), we get the following equation for \(q_1\):

$$\begin{aligned} q_1'(z) = -2\frac{q_0(z)}{K_0(z)}q_1(z) + \frac{q_0(z)^2 K_1(z)}{K_0(z)^2}, \quad z\in (-h,0), \end{aligned}$$
(3.6)

with “initial” condition

$$\begin{aligned} q_1(-h) = 0, \end{aligned}$$
(3.7)

in view of (3.1b) and (3.5b). The solution to (3.6)–(3.7) is

$$\begin{aligned} q_1(z) = \exp \left( -2\int \limits _{-h}^z\frac{q_0(y)}{K_0(y)}\,\mathrm {d}y\right) \int \limits _{-h}^z \exp \left( 2\int \limits _{-h}^s \frac{q_0(y)}{K_0(y)}\,\mathrm {d}y\right) \frac{q_0(s)^2 K_1(s)}{K_0(s)^2}\,\mathrm {d}s.\end{aligned}$$

Thus

$$\begin{aligned} q(0)&= q_0(0) + \varepsilon \int \limits _{-h}^0\exp \left( -2\int \limits _{s}^0 \frac{q_0(y)}{K_0(y)}\,\mathrm {d}y\right) \frac{q_0(s)^2 K_1(s)}{K_0(s)^2}\,\mathrm {d}s+ o(\varepsilon ) \\&= q_0(0)\left[ 1+\frac{\varepsilon }{q_0(0)}\int \limits _{-h}^0\exp \left( -2\int \limits _{s}^0 \frac{q_0(y)}{K_0(y)}\,\mathrm {d}y\right) \frac{q_0(s)^2 K_1(s)}{K_0(s)^2}\,\mathrm {d}s\right] + o(\varepsilon ). \end{aligned}$$

The formula

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}\varepsilon }\arg (1+\varepsilon \xi )\bigg \vert _{\varepsilon =0} = \mathrm {Im}(\xi ), \quad \xi \in \mathbb {C}\end{aligned}$$

(cf. [1, 2]) enables us to express the change of \(\arg (q(0))\) due to the perturbation as

$$\begin{aligned} \mathrm {Im}\left\{ \frac{1}{q_0(0)}\int \limits _{-h}^0\exp \left( -2\int \limits _{s}^0 \frac{q_0(y)}{K_0(y)}\,\mathrm {d}y\right) \frac{q_0(s)^2 K_1(s)}{K_0(s)^2}\,\mathrm {d}s\right\} =:C_{K_1}. \end{aligned}$$
(3.8)

Therefore, by (3.2), a positive/negative value of \(C_{K_1}\) corresponds to a decrease/increase of \(\arg (\psi (0))\) from the reference value \(\arg (\psi _0(0))\). Moreover, (3.3) implies

$$\begin{aligned} \exp \left( -2\int \limits _{s}^0 \frac{q_0(y)}{K_0(y)}\,\mathrm {d}y\right) = \frac{\psi _0(s)^2}{\psi _0(0)^2} \quad \text {and} \quad \frac{q_0(s)}{K_0(s)} = \frac{\psi _0'(s)}{\psi _0(s)},\end{aligned}$$

thus formula (3.8) can be rewritten in terms of \(\psi _0\) as

$$\begin{aligned} C_{K_1} = \frac{1}{K_0(0)}\,\mathrm {Im}\left\{ \frac{1}{\psi _0(0)}\int \limits _{-h}^0(\psi _0'(s))^2 K_1(s)\,\mathrm {d}s\right\} . \end{aligned}$$
(3.9)

4 Discussion

The easiest situation where formula (3.9) can be applied is the case of constant \(K_0 \equiv 1\) on \((-\infty ,0]\) (cf. the discussion in [1, 2]). Here we have

$$\begin{aligned} \psi _0(z) = \frac{1-\mathrm {i}}{2}\mathrm {e}^{(1+\mathrm {i})z},\end{aligned}$$

hence, according to (3.2) and (3.9), the change of \(\arg (\psi (0))\) is given by

$$\begin{aligned} -C_{K_1}&= -\mathrm {Im}\left\{ (1+\mathrm {i})\int \limits _{-h}^0\mathrm {e}^{2(1+\mathrm {i})s}K_1(s)\,\mathrm {d}s\right\} \nonumber \\&= -\sqrt{2}\int \limits _{-h}^0\mathrm {e}^{2s}\sin \left( 2s+\frac{\pi }{4}\right) K_1(s)\,\mathrm {d}s, \end{aligned}$$

which corresponds, up to a scaling, to the formula derived in [2] applied to an upper-layer perturbation which is confined to the interval \((-h,0)\).

It is interesting to notice that the presence of a sine function in the integrand is not just a feature of the special case of constant \(K_0\). In fact, if we write

$$\begin{aligned} \psi _0(0) = \vert \psi _0(0)\vert \mathrm {e}^{\mathrm {i}\arg (\psi _0(0))} \quad \text {and} \quad \psi _0'(s) = \vert \psi _0'(s)\vert \mathrm {e}^{\mathrm {i}\arg (\psi _0'(s))},\end{aligned}$$

then formula (3.9) becomes

$$\begin{aligned} C_{K_1}&= \frac{1}{K_0(0)\vert \psi _0(0)\vert }\mathrm {Im}\left\{ \int \limits _{-h}^0\vert \psi _0'(s)\vert ^2\mathrm {e}^{\mathrm {i}[2\arg (\psi _0'(s))-\arg (\psi _0(0))]}K_1(s)\,\mathrm {d}s\right\} \\&= \frac{1}{K_0(0)\vert \psi _0(0)\vert }\int \limits _{-h}^0\vert \psi _0'(s)\vert ^2\sin [2\arg (\psi _0'(s))-\arg (\psi _0(0))]K_1(s)\,\mathrm {d}s. \end{aligned}$$

From this it is apparent that in general it is very difficult to tell at first glance whether a given perturbation \(K_1\) of \(K_0\) will lead to a positive or negative increment of the surface deflection angle. For instance, suppose that \(K_1(z)>0\) for each \(z\in (-h,0)\). Due to the periodicity of the sine, we may assume that \(2\arg (\psi _0'(s))-\arg (\psi _0(0)) \in (-\pi ,\pi ]\) for each \(s\in (-h,0)\). Let us denote

$$\begin{aligned} P&:=\{s\in (-h,0) : 2\arg (\psi _0'(s))-\arg (\psi _0(0)) \in (0,\pi )\}, \\ N&:=\{s\in (-h,0) : 2\arg (\psi _0'(s))-\arg (\psi _0(0)) \in (-\pi ,0)\}. \end{aligned}$$

Since \(\sin [2\arg (\psi _0'(s))-\arg (\psi _0(0))]\) is positive for \(s\in P\) and negative for \(s\in N\), the sign of \(C_{K_1}\) will depend both on the measures \(\vert P\vert \) and \(\vert N\vert \) and on the relative sizes of the values attained by \(K_1\) on P and N.

Finally, let us look at one last example that shows how the “Riccati” variant (3.8) of (3.9) can be useful in its own way, not just as a step in the derivation of (3.9). This is especially true if the solution is found using the reformulation (3.1) in the first place. In [9] the solution to the Riccati formulation of the problem for the case

$$\begin{aligned} K_0(z) = [3(z+h)+1]^{4/3}, \quad z\in (-h,0)\end{aligned}$$

is determined to be

$$\begin{aligned}q_0(z) = -S(z)-(1-\mathrm {i})S(z)^2\tan [(1-\mathrm {i})S(z)+C],\end{aligned}$$

where \(S(z)=[3(z+h)+1]^{1/3}\) and

$$\begin{aligned} C = -1+\mathrm {i}+\frac{\mathrm {i}}{2}\log \left( \frac{1-\mathrm {i}}{\mathrm {i}-5}\right) .\end{aligned}$$

Using the identity

$$\begin{aligned} \tan (\theta ) = -\frac{\mathrm {d}}{\mathrm {d}\theta }\log (\cos (\theta ))\end{aligned}$$

and the fundamental theorem of calculus, a straightforward computation yields

$$\begin{aligned} \exp \left( -2\int \limits _s^0\frac{q_0(y)}{K_0(y)}\,\mathrm {d}y\right) = \left( \frac{S(0)\cos [(1-\mathrm {i})S(s)+C]}{S(s)\cos [(1-\mathrm {i})S(0)+C]}\right) ^2,\end{aligned}$$

which can be directly plugged into (3.8) without first having to determine and differentiate \(\psi _0\).

We conclude by saying that formula (3.9) may be of some relevance especially in view of the fact that, as pointed out in the introduction, analytic solutions are hard to find and the few ones available show an intricate behaviour even for relatively simple eddy viscosity profiles. The perturbation approach that we have pursued here may be useful to bring some new insight into this delicate problem, but it would be interesting to obtain further results that do not require the smallness assumption that is intrinsic to the perturbation method. On the other hand, it would also be of great importance to bring time dependence into the picture; our current knowledge of time-dependent Ekman flows is very limited and is restricted to some special cases (see [13]), but it seems likely that the technique used in this paper may be employed to obtain further results in this direction.