1 Introduction

The Earth’s atmosphere can be divided into several layers based on the behaviour of its temperature (see [19]). These layers are, starting from ground level upwards, the troposphere, the stratosphere, the mesosphere and the thermosphere. A further region, beginning about 500 km above the ground level, is the exosphere, which fades away into the realm of interplanetary space. The troposphere, whose height decreases poleward (from about 17–18 km at the Equator to about 7–8 km at the poles), contains more than 75% of all of the air in the atmosphere, and almost all of the water vapour (which forms clouds and rain). This is the region where the familiar weather phenomena occur. Throughout the troposphere there is a decrease in temperature with height, due to a decreasing pressure, since air is a compressible fluid. The lowest part—roughly the lower third—of the troposphere is called the atmospheric boundary layer, and it is here that friction plays an important role, while higher up, from the stratosphere upwards, the air flow is practically inviscid. For a better understanding of the flow dynamics, it is useful to divide the planetary boundary layer in two parts: the surface layer and the Ekman layer. Within the surface layer, confined to the lowest few meters of the atmosphere (above the ground), the velocity profile is adjusted so that the horizontal frictional stress is nearly independent of height. In contrast to this, in the Ekman layer, located on top of the surface layer and extending to a height of about 1 km, on average, the flow is governed by a three-way balance among the Coriolis force (due to the Earth’s rotation around its polar axis), the pressure gradient force, and the viscous stress (see [12, 23]). Primarily the air flow is horizontal (the horizontal velocities are typically of the order of 10 m s\(^{-1}\), about a factor \(10^4\) larger than the vertical velocity—see [24]), and in such a situation friction between the main horizontal motion and the bottom boundary is unavoidable. Friction reduces the velocity in the vicinity of the bottom, thus creating a vertical shear, with the main variations of the horizontal velocity in the vertical direction and controlled by the eddy viscosity, which captures the frictional effects on time scales less than a day. The Coriolis effect also plays a role, deflecting (in regions away from the Equator) the orientation of the flow, an aspect that differs greatly from the boundary-layer in non-rotating fluids (also, at the Equator the Coriolis force does not play a role—see the discussion of the ocean context in [1, 3, 5, 9, 11, 15, 17, 20]). The theory of the flow in the atmospheric boundary layer at mid-latitudes was pioneered in 1908 by Åkerblom, who adapted to the atmosphere an approach developed by Ekman in 1902 for motion in the upper layers of the ocean under the influence of a steady wind (see the discussion in [10]). For constant eddy viscosity there is an explicit solution, but this assumption is too restrictive. The dynamics of the atmospheric boundary-layer is very important in applications, for example, other than meteorology (weather prediction and climate studies), in the control and management of air pollution (since the dispersal of smog in urban environments depends strongly on meteorological conditions) and in agriculture (e.g. dewfall and frost formation). For this reason, it is important, both from the theoretical as well as from the practical point of view, to understand the flow dynamics of the atmospheric boundary-layer in the context of height-dependent eddy viscosities.

The explicit classical solution, valid for constant eddy viscosity, has no general counterpart for height-dependent eddy viscosities: only a few cases are known—see the discussion in [7, 13, 14, 18, 21]—and one typically performs numerical simulations. Note that while recently it was established that in general the Ekman flow pattern is spiralling and the horizontal speed is monotone with height (see [7]), there remains the important issue of the deflection angle. Apart from some special, explicit solutions, there is little to go by, whether in the context of atmospheric flows or regarding wind-generated ocean currents (even though for the latter, a perturbative approach, developed in [2, 4] and applicable to small perturbations of a constant eddy viscosity, can be used to predict whether the deflection angle is smaller or bigger than the 45\(^\circ \) of the classical solution, valid for constant eddy viscosity). For this reason, it is important to identify new classes of height-dependent eddy viscosities that permit explicit calculations.

In this paper we describe an approach that combines a Liouville transformation of the governing equations for atmospheric Ekman flows with a suitable substitution to reduce the problem to a Riccati ordinary differential equation. This way, starting with the Riccati equation, we can work backwards and identify classes of height-dependent eddy viscosities for which explicit calculations are possible. In Sect. 2 we recall the governing equations and in Sect. 3 we implement our approach of reduction to a Riccati equation. Section 4 is devoted to a simple illustration of the implementation of the method.

2 Preliminaries

The governing equations for mesoscale steady air flow at mid-latitudes in the Ekman layer are (see the discussion in [16])

$$\begin{aligned} f(u-u_g)= & {} \frac{\mathrm{d}}{\mathrm{d}z}\Big (K\,\frac{\mathrm{d}v}{\mathrm{d}z}\Big )\,, \end{aligned}$$
(2.1)
$$\begin{aligned} f(v-v_g)= & {} -\,\frac{\mathrm{d}}{\mathrm{d}z}\Big (K\,\frac{\mathrm{d}u}{\mathrm{d}z}\Big )\,, \end{aligned}$$
(2.2)

where (u, v) is the horizontal z-dependent wind velocity (z being the height), with zonal (West-to-East, in the sense of the Earth’s rotation) component u and meridional (positive meaning towards the North Pole) component v, \(u_g\) and \(v_g\) are the corresponding geostrophic wind components, while \(K=K(z)\) is the height-dependent eddy viscosity. Here

$$\begin{aligned} f=2\mathrm {\Omega }\sin \theta \end{aligned}$$

is the Coriolis parameter at fixed latitude, \(\theta \), in the Northern Hemisphere, while \(\mathrm {\Omega } \approx 7.29 \times 10^{-5}\) s\(^{-1}\) is the angular speed of rotation of the Earth and \(\theta \in (0,\pi /2]\) is the angle of latitude in right-handed rotating spherical coordinates (with \(\theta =0\) corresponding to the Equator and \(\theta =\pi /2\) to the North Pole). The boundary conditions for the system (2.1)–(2.2) are

$$\begin{aligned}&u=v=0 \quad \text {at}\quad z=0\,, \end{aligned}$$
(2.3)
$$\begin{aligned}&u \rightarrow u_g \quad \text {and}\quad v \rightarrow v_g\quad \text {for}\quad z \rightarrow \infty \,, \end{aligned}$$
(2.4)

expressing the fact that, due to the frictional properties of the flow below the Ekman layer, a no-slip condition holds at the bottom \(z=0\) of the layer, while at the top of the Ekman layer the horizontal components of the wind must be in geostrophic balance: above the Ekman layer the flow is geostrophic (pressure-driven). We note that the governing equations (2.1)–(2.2) are derived within the framework of the f-plane approximation, but one can easily adapt to the atmospheric setting the considerations made in [6, 8] for ocean flows to obtain a quite similar system in rotating spherical coordinates. Since it is standard procedure to work with the f-plane, we will consider in this paper the system (2.1)–(2.2) with the boundary conditions (2.3)–(2.4).

For constant eddy viscosity, the explicit classical solution is

$$\begin{aligned} \left\{ \begin{array}{l} u(z)=u_g - e^{-\gamma z}[u_g \cos (\gamma z)+v_g \sin (\gamma z)]\,,\\ v(z)=v_g + e^{-\gamma z}[u_g\sin (\gamma z) -v_g\cos (\gamma z)]\,, \end{array}\right. \end{aligned}$$
(2.5)

with

$$\begin{aligned} \gamma =\sqrt{\frac{f}{2K}}, \end{aligned}$$

showing the spiraling behaviour of the horizontal ageostrophic wind \((u-u_g,v-v_g)\), decreasing upwards: we can express the solution (2.5) in complex notation as

$$\begin{aligned}{}[u(z)-u_g]+i[v(z)-v_g]=- e^{-(1+i)\gamma z}\{u_g+iv_g\}\,, \end{aligned}$$

to visualize better the spiralling pattern diminishing with increasing height z. Of practical interest are not only constant eddy viscosities K, but continuous and bounded functions

$$\begin{aligned} K: [0,\infty ) \rightarrow [k_-,k_+]\,, \quad \text {where}\quad k_\pm ,\,k^*>0 \quad \text {are positive constants}\,, \end{aligned}$$

with

$$\begin{aligned} K(z) \rightarrow k^*\quad \text {as}\quad z \rightarrow \infty \quad \text {at a fast rate}, \end{aligned}$$

for some \(k^*\in [k_-,k_+]\), in the sense that there exists constants \(a,b,c>0\) such that

$$\begin{aligned} |K(z) - k^*| \le \frac{a}{1+bz^{2+c}}\,,\qquad z \ge 0\,. \end{aligned}$$

In particular, these restrictions can clearly allow profiles K(z) that are non-monotonic, these being especially relevant physically (see the discussion in [14]). As pointed out in [7], the Liouville substitution

$$\begin{aligned} s=\xi (z)= \int _0^z \frac{1}{K(s)}\,\mathrm{d}s\,,\qquad s \in [0,\infty )\,, \end{aligned}$$
(2.6)

transforms (2.1)–(2.2) to the second-order linear differential equation

$$\begin{aligned} \mathrm{\Psi }''(s) = i\,\alpha (s)\,\mathrm{\Psi }(s)\,,\qquad s>0\,, \end{aligned}$$
(2.7)

for the complex-valued function

$$\begin{aligned} \mathrm{\Psi }(s)=U(s)+iV(s) \end{aligned}$$
(2.8)

where we denoted

$$\begin{aligned} \left\{ \begin{array}{l} U(s)=u(z)-u_g\,,\\ V(s)=v(z)-v_g \,, \end{array}\right. \quad \text {for}\quad s=\xi (z)\,, \end{aligned}$$
(2.9)

with

$$\begin{aligned} \alpha (s)=f\,K(z) \quad \text {for}\quad s=\xi (z) \ge 0\,, \end{aligned}$$
(2.10)

and the prime denotes differentiation with respect to the s variable. The boundary conditions (2.3)–(2.4) are transformed into

$$\begin{aligned}&\mathrm{\Psi }=-(u_g + iv_g)\quad \text {at}\quad s=0\,, \end{aligned}$$
(2.11)
$$\begin{aligned}&\mathrm{\Psi } \rightarrow 0 \quad \text {for}\quad s \rightarrow \infty \,. \end{aligned}$$
(2.12)

One can easily verify that the function \(s \mapsto \alpha (s)\) captures the monotonicity properties of the function \(z \mapsto K(z)\). As showed in [7], in the Northern Hemisphere, for any positive continuous function \(\alpha : [0,\infty ) \rightarrow (0,\infty )\) such that

$$\begin{aligned} \int _0^\infty s |\alpha (s)-k^*f|\,\mathrm{d}s < \infty \,, \end{aligned}$$
(2.13)

there is a unique solution to the boundary-value problem (2.7)–(2.11)–(2.12), obtained iteratively; the only difference for the Southern Hemisphere is that therein the Coriolis parameter \(f=2\Omega \sin \theta <0\) because \(\theta \in [-\frac{\pi }{2},0)\), and instead of f we should write |f| in the corresponding condition (2.13).

With the existence of solutions guaranteed under the mild (and physically realistic) condition (2.13), the issue that we address in the next section is that of providing a method that yields explicit solutions to the second-order linear differential equation (2.7) with boundary conditions (2.11)–(2.12).

3 Main Result

It is convenient to perform a scaling of the Sturm–Liouville problem (2.7)–(2.11)–(2.12), to reduce the number of parameters. Setting

$$\begin{aligned} \mathrm{\Phi }(s)=-\,\frac{\Psi (s)}{u_g + iv_g}\,,\qquad s>0\,, \end{aligned}$$
(3.1)

we transform (2.7)–(2.11)–(2.12) into the equivalent Sturm–Liouville problem

$$\begin{aligned}&\mathrm{\Phi }''(s) = i\,\alpha (s)\,\mathrm{\Psi }(s)\,,\qquad s>0\,, \end{aligned}$$
(3.2)
$$\begin{aligned}&\mathrm{\Phi }=1 \quad \text {at}\quad s=0\,, \end{aligned}$$
(3.3)
$$\begin{aligned}&\mathrm{\Phi } \rightarrow 0 \quad \text {for}\quad s \rightarrow \infty \,. \end{aligned}$$
(3.4)

Let us now prove the following result.

Theorem 3.1

The solution \(\mathrm{\Phi }\) of the Sturm–Liouville problem (3.2)–(3.3)–(3.4) can be expressed in the form

$$\begin{aligned} \mathrm{\Phi }(s)=e^{ \int _0^s h(\tau )\,\mathrm{d}\tau }\,,\qquad s \ge 0\,, \end{aligned}$$
(3.5)

where h(s) is a solution of the Riccati equation

$$\begin{aligned} h'(s)+h^2(s)=i\alpha (s)\,,\qquad s>0\,, \end{aligned}$$
(3.6)

with

$$\begin{aligned} \lim \limits _{s \rightarrow \infty } h(s)=-\,(1+i)\sqrt{\frac{k^*f}{2}}\,. \end{aligned}$$
(3.7)

Proof

Let us first show that \(\mathrm{\Phi }(s)>0\) for all \(s \ge 0\). If this were not so, sine \(\mathrm{\Phi }(0)=1\), by continuity there would exist \(s_0=\inf \{s>0:\ \mathrm{\Phi }(s)=0\}>0\) with \(\mathrm{\Phi }(s_0)=0\). But then, multiplying (3.2) by the complex conjugate \(\overline{\mathrm{\Phi }}(s)\) of \(\mathrm{\Phi }(s)\), an integration on \([s_0,\infty ]\) yields the contradiction

$$\begin{aligned} -\int _{s_0}^\infty |\mathrm{\Phi }'(\tau )|^2 \mathrm{d}\tau =i \int _{s_0}^\infty \alpha (\tau ) |\mathrm{\Phi }(\tau )|^2 \mathrm{d}\tau \,, \end{aligned}$$

since the right side is purely imaginary and the left side is real. To justify the above integral relation, note that the iterative process yielding \(\mathrm{\Psi }(s)=-[u_g+iv_g] \mathrm{\Phi }(s)\) ensures (see the discussion in [7])

$$\begin{aligned} \left\{ \begin{array}{l} \mathrm{\Psi }(s) \sim e^{-(1+i)s \sqrt{k^*f/2}}\\ \mathrm{\Psi }'(s) \sim -\,(1+i)\sqrt{\frac{k^*f}{2}}\,e^{-(1+i)s \sqrt{k^*f/2}} \end{array}\right. \qquad \text {for}\qquad s \rightarrow \infty \,. \end{aligned}$$

In combination with (2.13), the above asymptotic behaviour validates the integral relation. As a byproduct, we see that \(\lim \limits _{s \rightarrow \infty }\mathrm{\Phi }(s)=0\) with

$$\begin{aligned} \lim \limits _{s \rightarrow \infty } \frac{\mathrm{\Phi '}(s)}{\mathrm{\Phi }(s)}=-\,(1+i)\sqrt{\frac{k^*f}{2}}\,. \end{aligned}$$
(3.8)

The above considerations prove that \(\mathrm{\Phi }(s)>0\) for all \(s \ge 0\). This permits us to introduce the function

$$\begin{aligned} h(s)=\frac{\mathrm{\Phi }'(s)}{\mathrm{\Phi }(s)}\,,\qquad s \ge 0\,. \end{aligned}$$

Using (3.2), it is easy to check that h solves (3.6). Taking now (3.3) into account, an integration of the relation displayed above leads us to the formula (3.5). Finally, the asymptotic behaviour (3.7) is simply another way to write (3.8).

4 Example

While a general solution of the Riccati equation (3.6) is not available explicitely, there are various types of cases in which the solutions are available in closed form. This fact can be used to construct explicit solutions to the atmospheric flow boundary-value problem (3.2)–(3.3)–(3.4). We now illustrate the procedure by means of an example for which we present the full details. But let us first present some general observations. Note that (2.13) does not accommodate polynomial expressions in s for the function \(\alpha (s)\), so that to the Airy and Hermite transcendent functions (corresponding to polynomials of degree 1 and 2, respectively, in the Riccati equation (3.6)) are not relevant. However, a general property of (3.6) is useful, namely the fact that one available particular solution allows to construct its general solution. Indeed, if \(h_1\) is the available solution, rewriting (3.6) on an interval where \(hh_1 \ne 0\) as

$$\begin{aligned} \frac{h'}{h^2}=-1 + i\,\frac{\alpha }{h}\,, \end{aligned}$$

the change of variables \(y=1/h\) leads us to the equivalent problem

$$\begin{aligned} y' + i \alpha y -1=0\,. \end{aligned}$$

Since \(y_1=1/h_1\) is a solution of the above first-order non-homogeneous linear ordinary differential equation, we see that \(Y=y-y_1\) solves

$$\begin{aligned} Y' + i \alpha Y =0\,, \end{aligned}$$

so that

$$\begin{aligned} Y(s)=Y_0\, e^{-i \int _0^s \alpha (\tau )\,\mathrm{d}\tau }\,,\qquad s \ge 0\,, \end{aligned}$$

for some initial data\(Y_0 \in {\mathbb C}\). Consequently \(y=Y+y_1\), and thus

$$\begin{aligned} h(s)=\frac{h_1(s)}{1 + h_1(s)Y(s)}=\frac{h_1(s)}{1 + h_1(s)Y_0\, e^{-i \int _0^s \alpha (\tau )\,\mathrm{d}\tau }}\,. \end{aligned}$$

Let us now discuss a detailed class of examples. We claim that if \(\beta : [0,\infty ) \rightarrow {\mathbb R}\) is a twice-continuously differentiable function that approaches a constant at infinity and satisfies

$$\begin{aligned} \beta '(s) + \beta ^2(s) \ge 0\,,\qquad s \ge 0\,, \end{aligned}$$

(e.g. \(\beta (s)=d \tanh (\lambda s)\) with \(d \ge \lambda >0\) is of this type), then

$$\begin{aligned} h(s)=\beta (s) + i \sqrt{\beta '(s) + \beta ^2(s)}\,,\qquad s \ge 0\,, \end{aligned}$$

solves (3.6) for

$$\begin{aligned} \alpha (s)=\frac{\beta ''(s)+ 2\beta (s)\beta '(s)}{2\sqrt{\beta '(s) + \beta ^2(s)}} + 2\beta (s) \sqrt{\beta '(s) + \beta ^2(s)}\,,\qquad s \ge 0\,. \end{aligned}$$

While this claim can be verified directly, an alternative approach is to note that if \(h=\beta +i \eta \) is a solution of (3.6), with \(\beta \) and \(\eta \) real-valued functions, then \(\alpha =\eta ' +2\beta \eta \) (being a real-valued function) and \(\beta '+\beta ^2-\eta ^2=0\). The latter equation determines uniquely the real-valued function \(\eta \), and leads us to the above formulas.

For more detailed information on explicit solution of the Riccati equations we refer to [22].