2.1. Inviscid solution

The base flow $\boldsymbol {u}^{\left (0\right )}$ forced by precession can be found by solving the Navier–Stokes equations (2.1a,b) at first order in $\varepsilon$. Following Lagrange et al. (Reference Lagrange, Meunier and Eloy2016), we only keep the first term on the right-hand side of (2.1). The solution $\boldsymbol {u}^{\left (0\right )}$ is then composed of two parts: a particular solution in the form of a horizontal shear flow and a sum of Kelvin modes of azimuthal wavenumber $m=1$, which are the solutions of the homogeneous equation. Thus

(2.2)\begin{equation} \boldsymbol{u}^{\left(0\right)} = \varepsilon \left( \boldsymbol{u}_{{{shear}}} + \sum_{j = 1}^{\infty} A_j {\textbf{v}} _j \right) \exp({\text{i}(\omega t + m\theta)}) + \text{c.c.}, \end{equation}

with $\text {c.c.}$ standing for ‘complex conjugate’. In (2.2), $\boldsymbol {u}_{{{shear}}} = \omega z(\text {i} \hat {\boldsymbol {r}} -\hat {\boldsymbol {\theta }})/(2 - \omega )$ is the horizontal shear, $A_j$ are the amplitudes of Kelvin modes of axial wavenumber $k_j = (2j-1){\rm \pi} /h$, and ${\textbf {v}} _j$ their associated velocity fields given by

(2.3) \begin{equation} {\textbf{v}} _j = \left(\begin{array}{@{}r@{}} \text{i} u \sin \left( k_j z \right) \\ v \sin \left( k_j z \right) \\ \text{i} w \cos \left( k_j z \right)\end{array}\right), \end{equation}

with $u$, $v$ and $w$, real functions of $r$,

(2.4a)$$\begin{gather} u(r) = \frac{ - 1}{4 - \omega ^{2} }\left( \omega \delta J'_m \left(\delta r\right) + \frac{2m}{r}J_m \left( \delta r \right) \right), \end{gather}$$

(2.4b)$$\begin{gather}v(r) = \frac{1}{4 - \omega^{2}}\left( 2\delta J'_m \left( \delta r \right) + \frac{\omega m}{r}J_m \left( \delta r \right) \right), \end{gather}$$

(2.4c)$$\begin{gather}w(r) = \frac{k}{\omega}J_m \left(\delta r\right), \end{gather}$$

with $J_m$ the Bessel function of the first kind, $\delta = k_j(4/\omega ^{2} -1)^{1/2}$ playing the role of a radial wavenumber, and primes noting differentiation.

Here, we will assume that the forcing is resonant. This occurs when the forcing frequency takes particular values, $\omega =\omega _{n,j}$, such that $u(1)=0$ for a particular Kelvin mode. In that case, this Kelvin mode amplitude $A_j$ is asymptotically large and the radial wavenumber takes the value $\delta = \delta _{n,j} = k_j(4/\omega _{n,j}^{2} -1)^{1/2}$. Each resonance is associated with a couple of strictly positive integers $(n,j)$: $n$ numbers the radial wavenumber $\delta$ (the function $u(r)$ has exactly $n$ zeros for $0< r\le 1$), and $j$ numbers the axial wavenumber $k_j$ such that $k_j = (2j-1){\rm \pi} /h$. In Meunier et al. (Reference Meunier, Eloy, Lagrange and Nadal2008), we also referred to these resonances by the phrase ‘$j$th resonance of the $n$th Kelvin mode’.

At resonance, the amplitude $A_j$ of the resonant Kelvin mode can be calculated by either invoking viscous effects, in which case $A_j=O({ {Ek}}^{-1/2})$, or by invoking nonlinear effects, in which case $A_j=O(\varepsilon ^{-2/3})$ (Meunier et al. Reference Meunier, Eloy, Lagrange and Nadal2008). In either case, the base flow is dominated by this single resonant Kelvin mode of amplitude $A=A_j$,

(2.5)\begin{equation} \boldsymbol{u}^{\left(0\right)} \approx A {\textbf{v}} \exp({\text{i}(\omega t + m \theta)}) + \text{c.c.} , \end{equation}

whose velocity field ${\textbf {v}}$ is given by (2.3) and we have dropped the indices $j$ for simplicity.

2.2. Viscous solution

The base flow given in (2.5) is a solution of the inviscid problem. To satisfy the viscous boundary condition on the walls, it has to be complemented by a boundary layer solution confined in the wall regions of thickness $O({ {Ek}}^{1/2})$. The total flow at leading order is then

(2.6)\begin{equation} \boldsymbol{u}^{(0)}_{{tot.}} = A\left({\textbf{v}} + \widetilde{{\textbf{v}}} \right) \exp({\text{i}(\omega t + m \theta)}) +\text{c.c.} \end{equation}

The viscous solution $\widetilde {{\textbf {v}}}$ can be expressed as a series in powers of ${ {Ek}}^{1/2}$. This viscous solution has a non-zero axial component $\tilde w$ on the end walls at order ${ {Ek}}^{1/2}$, forced through Ekman pumping (Meunier et al. Reference Meunier, Eloy, Lagrange and Nadal2008),

(2.7)\begin{equation} \tilde w (z={-}h/2) = {{Ek}}^{1/2} k^{2} \frac{\text{i}-1}{2\sqrt{2}\omega^{2}} \left(\frac{2+\omega }{(2 - \omega)^{1/2}} +\text{i}\frac{2-\omega}{(2 + \omega)^{1/2}}\right) J_m (\delta r) \sin(kh/2), \end{equation}

with an opposite flow on the upper wall in $z=h/2$.

To balance these non-zero axial flows on the end walls, the inviscid solution has to be corrected at order ${ {Ek}}^{1/2}$. A simple way to achieve this correction is to write the axial wavenumber as a series in powers of ${ {Ek}}^{1/2}$:

(2.8)\begin{equation} k = k^{(0)} + {{Ek}}^{1/2} k^{(1)} + \cdots,\quad \mbox{with } k^{(0)} = (2j-1){\rm \pi}/h \end{equation}

and

(2.9)\begin{equation} k^{(1)} = \frac{k^{(0)}}{\sqrt{2}\omega h} (1 + \text{i}) \left(\frac{2+\omega }{(2 - \omega)^{1/2}} + \text{i}\frac{2-\omega}{(2 + \omega)^{1/2}}\right). \end{equation}

A similar correction exists for the radial wavenumber $\delta$ since it is proportional to $k$. When there is no ambiguity, we will refer to $k^{(0)}$ as $k$ in the following.

4.1. Interaction between equatorial Coriolis force and the Kelvin mode

The first type of flow that forces a geostrophic flow is ${\textbf {v}}_k$, the axisymmetric and steady solution of the linearised Navier–Stokes equation forced by the equatorial Coriolis force, i.e. the second term on the right-hand side of (2.1). It thus satisfies the following equations:

(4.1a,b)\begin{equation} 2\hat{\boldsymbol{z}}\times {\textbf{v}}_{k} + \boldsymbol{\nabla} p_{k} = \boldsymbol{F}_{k},\quad \boldsymbol{\nabla}\cdot{\textbf{v}}_{k} = 0,\quad\mbox{with } \boldsymbol{F}_{k} = \varepsilon A \,\boldsymbol{d} \times {\textbf{v}} \end{equation}

and $\boldsymbol {d} = -\hat {\boldsymbol {r}} + \text {i} \hat {\boldsymbol {\theta }}$. A particular solution of this forced equation gives an axial component of velocity

(4.2)\begin{equation} w_k ={-}\text{i} \varepsilon A \frac{1}{2kr}\frac{\text{d}(r w)}{\text{d} r} \sin (kz). \end{equation}

In order to cancel the resulting inflow on the end walls, a geostrophic flow is necessary such that its viscous part satisfies $\tilde {w}_0 = -w_k$ on the lower wall. Using (3.7), it then leads to a geostrophic flow of order ${ {Ek}}^{-1/2}$ higher than ${\textbf {v}}_k$ itself, with azimuthal velocity

(4.3)\begin{equation} \left(\boldsymbol{u}_0\right)_{k} ={-}\varepsilon {{Ek}}^{{-}1/2} \mbox{Im}(A) \frac{2 w}{k} \sin (kh/2) \hat{\boldsymbol{\theta}}, \end{equation}

with $\sin (kh/2) = (-1)^{j-1}$, because $k=(2j-1){\rm \pi} /h$. Note that velocities noted with the bold letter ‘$u$’, such as $(\boldsymbol {u}_0)_{k}$, refer to the real component of the flow and include the complex conjugate, contrarily to velocities denoted by a bold ‘$v$’, such as ${\textbf {v}}_{k}$, for instance.

4.2. Interaction between the Kelvin mode and the shear flow

The second flow that forces a geostrophic flow is ${\textbf {v}}_{shear}$, the flow forced by the nonlinear interaction of the Kelvin mode with the shear flow. This flow is a solution of

(4.4a,b)\begin{equation} 2\hat{\boldsymbol{z}}\times {\textbf{v}}_{shear} + \boldsymbol{\nabla} p_{shear} = \boldsymbol{F}_{shear}, \quad \boldsymbol{\nabla}\cdot{\textbf{v}}_{shear} = 0, \end{equation}

with

(4.5)\begin{equation} \boldsymbol{F}_{shear} ={-}\varepsilon (A\bar{\boldsymbol{u}}_{shear} \cdot \boldsymbol{\nabla} {\textbf{v}} + \bar{A}\bar{{\textbf{v}}} \cdot \boldsymbol{\nabla} {\boldsymbol{u}}_{shear}). \end{equation}

The component $\boldsymbol {F}_{shear} \cdot \hat {\boldsymbol {\theta }}$ has a non-zero real part, which yields through the incompressibility condition to a non-zero axial component of ${\textbf {v}}_{shear}$,

(4.6)\begin{equation} w_{shear} = \varepsilon \, \mbox{Im}(A) \frac{1}{2k r}\frac{\partial(r q_{shear})}{\partial r}, \end{equation}

with

(4.7)\begin{equation} q_{shear} ={-} w(r) \sin(kz) + w(r) \frac{2 \delta \omega z}{(2-\omega)\sqrt{4-\omega^{2}}} \cos(kz) . \end{equation}

At the end walls $z=\pm h/2$, this axial flow balances exactly the axial flow $w_k$ given by (4.2). It means that a geostrophic flow is created that cancels $(\boldsymbol {u}_0)_{k}$,

(4.8)\begin{equation} \left(\boldsymbol{u}_0\right)_{shear} ={-}\left(\boldsymbol{u}_0\right)_{k}, \end{equation}

with $(\boldsymbol {u}_0)_{k}$ given by (4.3).

4.3. Interaction of the Kelvin mode with itself in the bulk

The third flow that forces a geostrophic flow is the mode ${\textbf {v}}_{2k}$, the axisymmetric and steady solution of the linearised Navier–Stokes equation forced by the nonlinear interactions of the Kelvin mode with itself, which satisfies

(4.9a,b)\begin{equation} 2\hat{\boldsymbol{z}}\times {\textbf{v}}_{2k} + \boldsymbol{\nabla} p_{2k} = \boldsymbol{F}_{2k},\quad \boldsymbol{\nabla}\cdot{\textbf{v}}_{2k} = 0,\quad\mbox{with } \boldsymbol{F}_{2k} ={-}|A|^{2} \bar{{\textbf{v}}} \cdot \boldsymbol{\nabla} {\textbf{v}}. \end{equation}

Note that the radial and axial components of the forcing are much larger than the azimuthal component, which vanishes for an inviscid Kelvin mode. The azimuthal component $\boldsymbol {F}_{2k} \cdot \hat {\boldsymbol {\theta }}$ is $O( |A|^{2} { {Ek}}^{1/2} )$ because of the viscous correction of the Kelvin mode given in (2.8). This azimuthal component of the forcing is treated, below, in § 4.3.2 through the component $(\boldsymbol {u}_0)_{Ek}$ of the geostrophic flow. We will first examine the geostrophic flow due to the radial and axial forcing.

4.3.1. Forcing at order $|A|^{2}$

At leading order $|A|^{2}$, a particular solution of (4.9a,b) has been given by Waleffe (Reference Waleffe1989),

(4.10)\begin{equation} {\textbf{v}}_{2k} ={-} |A|^{2} v_{2k}(r) \cos(kh) \cos\left( 2 k z\right)\hat{\boldsymbol{\theta}}, \quad \mbox{with } v_{2k}(r) = \frac{k}{\omega} \left(\frac{1}{2k}\frac{\text{d} (uv)}{\text{d} r} - vw \right), \end{equation}

with $u$, $v$ and $w$ the radial functions of the Kelvin mode, given in (2.4).

The velocity ${\textbf {v}}_{2k}$ is a solution of the inviscid problem. Although this solution seems to violate the statement of Greenspan (Reference Greenspan1969) that no geostrophic flow can be generated by an inviscid inertial mode, in fact it does not, since ${\textbf {v}}_{2k}$ is a non-geostrophic zonal flow (because of its dependence on $z$). The associated boundary layer solution can be found by the same method as the one described above for the geostrophic flow in § 3. It yields an inflow on the lower wall that can be written as

(4.11)\begin{equation} \tilde{w}_{2k}( z={-}h/2) = |A|^{2}{{Ek}}^{1/2} \frac{1}{2r} \frac{\text{d}\left( r v_{2k} \right)}{\text{d} r}. \end{equation}

To cancel this inflow, a geostrophic flow is necessary that will produce the opposite flow. This geostrophic flow is found using (3.7) and can be written as

(4.12)\begin{equation} \left(\boldsymbol{u}_0\right)_{2k} = 2|A|^{2} v_{2k}(r)\hat{\boldsymbol{\theta}} . \end{equation}

Note that the sum of the flow ${\textbf {v}}_{2k}$ and the above geostrophic flow gives a real component of the velocity

(4.13)\begin{equation} \boldsymbol{u}_{2k} + \left(\boldsymbol{u}_0\right)_{2k}= 2|A|^{2} v_{2k}(r) \left(1 - \cos(kh) \cos\left( 2 k z\right) \right)\hat{\boldsymbol{\theta}}, \end{equation}

that cancels for $z=\pm h/2$.

4.3.2. Forcing at order $|A|^{2} { {Ek}}^{1/2}$

At order $|A|^{2} { {Ek}}^{1/2}$, the inviscid Kelvin mode interacts nonlinearly with its viscous correction, which can be obtained by a series expansion of $k$ in powers of ${ {Ek}}^{1/2}$ given in (2.9). This nonlinear interaction can also be viewed as the interaction between the Kelvin mode and the flow resulting from Ekman pumping at the end walls. It yields a non-zero azimuthal forcing $\boldsymbol {F}_{2k} \cdot \hat {\boldsymbol {\theta }} = O( |A|^{2} { {Ek}}^{1/2} )$.

Through (4.9a,b), this forcing gives rise to a radial velocity $u_{Ek} = \frac {1}{2}\boldsymbol {F}_{2k} \cdot \hat {\boldsymbol {\theta }}$. This radial component induces an axial flow $w_{Ek}$ through the incompressibility condition

(4.14)\begin{equation} w_{Ek} (z={-}h/2) = |A|^{2} {{Ek}}^{1/2} \frac{1}{2r} \frac{\text{d}\left( r v_{Ek} \right)}{\text{d} r},\quad \mbox{with } v_{Ek}(r) = \frac{h }{\omega} u w\,\mbox{Im}\left(k^{(1)}\right), \end{equation}

where $k^{(1)}$ is the correction to the axial wavenumber given in (2.9) and $v_{Ek}$ is found by integrating $\boldsymbol {F}_{2k} \cdot \hat {\boldsymbol {\theta }}$ between $z=-h/2$ and 0 (see also (5.3), below).

To cancel this flow, a geostrophic flow $(\boldsymbol {u}_0)_{Ek}$ is necessary

(4.15)\begin{equation} \left(\boldsymbol{u}_0\right)_{Ek} = 2 |A|^{2} v_{Ek}(r)\hat{\boldsymbol{\theta}}. \end{equation}

4.4. Interactions of the Kelvin mode with itself in the boundary layers

Last, we consider the flow $\widetilde {{\textbf {v}}}_{BL}$ forced in the end wall boundary layers through nonlinear interactions of the viscous Kelvin mode with itself. Without loss of generality, we consider the lower wall boundary layer of thickness $O({ {Ek}}^{1/2})$. In this boundary layer, at leading order, the total flow ${\textbf {v}}^{(0)}_{{tot.}}$ can be written as the sum of the inviscid solution and its viscous correction

(4.16) \begin{equation} {\textbf{v}}^{(0)}_{{tot.}} = A \sin \left({-}kh/2\right) \left(\begin{array}{@{}c@{}} \text{i} u_{{tot.}} \\ v_{{tot.}} \\ {{Ek}}^{1/2} \text{i}\, w_{{tot.}} \end{array}\right) \exp({\text{i}(\omega t + m\theta)}), \end{equation}

where $u_{{tot.}}$, $v_{{tot.}}$ and $w_{{tot.}}$ are functions of $r$ and $\tilde {z}$, the rescaled vertical coordinate given by (3.3). These functions are given by (A1a–c) in Appendix A.

Through nonlinear interaction with itself, this flow acts as a steady and axisymmetric forcing of the linearised Navier–Stokes equations

(4.17a)$$\begin{gather} {{Ek}}^{{-}1/2} \frac{\partial \tilde{p}_{BL} }{\partial \tilde{z}} \hat{\boldsymbol{z}} + 2\hat{\boldsymbol{z}}\times \widetilde{{\textbf{v}}}_{BL} + \frac{\partial \tilde{p}_{BL} }{\partial r} \hat{\boldsymbol{r}} - \frac{\partial^{2} \widetilde{{\textbf{v}}}_{BL} }{\partial \tilde{z}^{2}} = \boldsymbol{F}_{BL}, \end{gather}$$

(4.17b)$$\begin{gather}{{Ek}}^{{-}1/2} \frac{\partial \tilde{w}_{BL} }{\partial \tilde{z}} + \frac{1}{r} \frac{\partial\left( r \tilde{u}_{BL} \right)}{\partial r} = 0, \end{gather}$$

where the forcing term can be written as

(4.18) \begin{equation} \boldsymbol{F}_{BL} (r,\tilde{z}) ={-}{\bar{\textbf{v}}}^{(0)}_{{tot.}} \cdot \boldsymbol{\nabla} {\textbf{v}}^{(0)}_{{tot.}} = |A|^{2} \sum_i \left( \begin{array}{@{}c@{}} a_i + c_i \tilde{z} \\ \text{i} b_i + \text{i} d_i \tilde{z} \\ O({{Ek}}^{1/2}) \end{array}\right) \text{e}^{-\kappa_i \tilde{z}}, \end{equation}

with the scalars $\kappa _i$ and the functions of $r$, $a_i$, $b_i$, $c_i$ and $d_i$ are given in Appendix A.

Equations (4.17a,b) can be solved by expanding $\widetilde {{\textbf {v}}}_{BL}$ in a series of powers of ${ {Ek}}^{1/2}$ and using the boundary conditions, $\widetilde {{\textbf {v}}}_{BL} = 0$ for $\tilde {z} =0$ and $\tilde {z}\to \infty$. The solution at relevant order is given by (A7) in Appendix A. This solution implies a non-zero axial flow $\tilde {w}_{BL}$ at order $|A|^{2} { {Ek}}^{1/2}$ on the lower end wall. To cancel this inflow, a geostrophic flow is necessary at order $|A|^{2}$ such that

(4.19)\begin{equation} \left(\boldsymbol{u}_0\right)_{BL} = 2 |A|^{2} v_{BL} (r)\hat{\boldsymbol{\theta}}, \end{equation}

with

(4.20)\begin{equation} v_{BL} (r) = \mbox{Re}\left[\sum_i \frac{a_i \kappa_i+\text{i} b_i (\kappa_i+2))}{\kappa_i (\kappa_i (\kappa_i+2)+2)} + 2 \frac{c_i (\kappa_i+1) \kappa_i^{2}+\text{i} d_i (\kappa_i (\kappa_i+2)^{2}+2)}{\kappa_i^{2} (\kappa_i (\kappa_i+2)+2)^{2}} \right] . \end{equation}

4.5. Total geostrophic flow

We now have calculated the five components of the geostrophic flows originating from five different inflows at the end walls. These components of the geostrophic flow, noted $(\boldsymbol {u}_0)_{k}$, $(\boldsymbol {u}_0)_{shear}$, $(\boldsymbol {u}_0)_{2k}$, $(\boldsymbol {u}_0)_{Ek}$ and $(\boldsymbol {u}_0)_{BL}$, are given, respectively, by (4.3), (4.8), (4.12), (4.15) and (4.19). Summing them gives the total forced geostrophic flow

(4.21)\begin{equation} \boldsymbol{u}_0 = \left(\boldsymbol{u}_0\right)_k + \left(\boldsymbol{u}_0\right)_{shear} + \left(\boldsymbol{u}_0\right)_{2k} + \left(\boldsymbol{u}_0\right)_{Ek} + \left(\boldsymbol{u}_0\right)_{BL} = A_0 v_0(r) \hat{\boldsymbol{\theta}}, \end{equation}

with

(4.22a,b)\begin{equation} A_0 = |A|^{2}, \quad v_0(r) = 2v_{2k}(r) + 2 v_{Ek}(r) + 2 v_{BL} (r), \end{equation}

and $v_{2k}$, $v_{Ek}$ and $v_{BL}$ given by (4.10), (4.14) and (4.20).

Note that the components $(\boldsymbol {u}_0)_k$ and $(\boldsymbol {u}_0)_{shear}$ exactly cancel each other, such that the amplitude of the geostrophic flow is proportional to $|A|^{2}$ and does not include any terms proportional to $\varepsilon A { {Ek}}^{-1/2}$. Because these two terms balance exactly, the resulting geostrophic flow does not depend on the precessional forcing anymore. The present results are thus applicable to any Kelvin mode, independently of the way it has been forced. In particular, as we shall see below, our calculations of the geostrophic flow is valid for Kelvin modes of arbitrary azimuthal wavenumber $m$.

4.6. Weakly nonlinear amplitude equations

So far, we have calculated the amplitude of the geostrophic flow $A_0$ in the steady regime, and we have found $A_0=|A|^{2}$. Based on this calculation and following the weakly nonlinear calculation in Meunier et al. (Reference Meunier, Eloy, Lagrange and Nadal2008), we can easily obtain the weakly nonlinear amplitude equations describing transient regimes.

To obtain the dynamic equations for $A$ and $A_0$, we have to correct the term describing nonlinear interaction between the geostrophic flow and the Kelvin mode. In Meunier et al. (Reference Meunier, Eloy, Lagrange and Nadal2008), we did not take into account the components $(\boldsymbol {u}_0)_{2k}$ and $(\boldsymbol {u}_0)_{Ek}$ of the geostrophic flow and there was also an error in the calculation of $(\boldsymbol {u}_0)_{Ek}$ (in particular, the term $w_4 \tilde z$ appearing in (A1c) was missing). At resonance, the dynamic amplitude equations are then given by

(4.23a)$$\begin{gather} \frac{\text{d} A}{\text{d} t} = \varepsilon \text{i} f - \left({{Ek}}^{1/2}\mu + {{Ek}}\nu\right) A + \text{i} \left(\sigma |A|^{2} + \xi A_0 \right) A, \end{gather}$$

(4.24b)$$\begin{gather}\frac{\text{d} A_0}{\text{d} t} = {{Ek}}^{1/2}\frac{2}{h}\left(|A|^{2} - A_0\right), \end{gather}$$

where the first four coefficients, $f$, $\mu$, $\nu$ and $\sigma$, can be calculated from Meunier et al. (Reference Meunier, Eloy, Lagrange and Nadal2008) and are given in table 1 for different resonances. The coefficient $f$ corresponds to the forcing by precession while $\mu$ and $\nu$ correspond to the Ekman and volume damping. The coefficient $\sigma$ corresponds to the nonlinear coupling of the Kelvin mode with the non-geostrophic flow ${\textbf {v}}_{2k}$ and with an unsteady elliptic flow noted ${\textbf {v}}_{2\omega }$ in Meunier et al. (Reference Meunier, Eloy, Lagrange and Nadal2008).

Table 1.

Coefficients appearing in the weakly nonlinear amplitude equations (4.23) for different resonances $(n,j)$ and different aspect ratios $h$.

The weakly nonlinear coefficient $\xi$, accounting for the nonlinear coupling between the Kelvin mode and the geostrophic flow, can be calculated as follows. We first calculate the nonlinear interaction of the Kelvin mode with the geostrophic flow: $\boldsymbol {u}_0 \boldsymbol {\cdot }\boldsymbol {\nabla } {\textbf {v}} + {\textbf {v}}\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}_0$. We then project this forcing onto the Kelvin mode using the natural Hermitian product over the cylinder volume. The coefficient $\xi$ is simply this Hermitian product normalised by the norm of the Kelvin mode (which is nothing else than its kinetic energy). It can be written $\xi = X/E$ with

(4.24)\begin{equation} X = 2 {\rm \pi}h |A|^{2} \int_0^{1} \left(\frac{2k}{\omega} u w v_0 - u v \frac{1}{r} \frac{\text{d} (r v_0)}{\text{d} r}\right)r\, \text{d} r, \end{equation}

with $E$ the kinetic energy of the mode

(4.25)\begin{equation} E = 2 {\rm \pi}h |A|^{2} \frac{2 k^{2} + m (2m - \omega)} {(4 - \omega^{2})\omega^{2}} J_m(\delta)^{2}, \end{equation}

and $v_0(r)$ given by (4.22a,b). The values of $\xi$ for the first nine resonances are plotted in figure 2. It shows that $\xi$ scales as $-h^{-4}$ for $h\ll 1$, changes sign for $h=O(1)$ and scales as $h$ for $h\gg 1$.

Figure 2.

Value of the weakly nonlinear coefficient $\xi$ as a function of the aspect ratio $h$. The value is given for different resonances $(n,j)$ as noted in the legend and we use the convention of a solid line for $\xi$ negative and a dashed line for $\xi$ positive.

This weakly nonlinear analysis will be used below to predict the amplitude $A_0$ of the geostrophic flow when nonlinear effects are not negligible.

5.1. Spatial structure of the forcing

In the Navier–Stokes equations, the interaction of the forced Kelvin mode with the time-dependent precession forcing gives rise to a forcing term $\boldsymbol {F}_k$ given in (4.1a,b). This term is noted $\mathcal {C}^{{inv}}$ in the paper of Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021). The three components of $\boldsymbol {F}_k$ along $\hat {\boldsymbol {r}}$, $\hat {\boldsymbol {\theta }}$ and $\hat {\boldsymbol {z}}$ are plotted in figure 3(a–c). They are in excellent agreement with the numerical results of Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) (see their figure 6d–f). The largest component is the $\theta$-component, which has a unique negative lobe with a minimum value equal to $-5.4 \times 10^{-5}$, while in the numerical simulations this value is approximately $-6 \times 10^{-5}$ (Albrecht et al. Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021). A small difference between the numerics and the theory is also found for the axial component: the numerical results have a 20 % smaller amplitude; they exhibit small undulations of the isocontours (probably due to inertial waves of small amplitudes emitted from the corners); in the simulations, this term vanishes at the top and the bottom because of the viscous boundary layers.

Figure 3.

Forcing terms of the different components of the geostrophic flow. The columns are, from left to right, the radial (a,d), azimuthal (b,e) and axial (c,f) components of the forcing terms. The components of the forcing term $\boldsymbol {F}_k$ given in (4.1a,b) are shown in panels (a–c). The components of the forcing term $\boldsymbol {F}_{\text {NL}} = \boldsymbol {F}_{shear} + \boldsymbol {F}_{2k} + \boldsymbol {F}_{BL}$ given in (4.5), (4.9a,b) and (4.18) are shown in panels (d–f). The dashed line shows the level line for zero. The parameters are $h=1.62$, ${ {Ek}}=1.5\times 10^{-4}$, $\omega =1.181$ and $\varepsilon =-0.00127$. To ease the comparison with figure 6 of Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021), we use the same contours.

The nonlinear forcing terms can be gathered into a single forcing term $\boldsymbol {F}_{\text {NL}} = \boldsymbol {F}_{shear} + \boldsymbol {F}_{2k} + \boldsymbol {F}_{BL}$, which is noted $\mathcal {R}^{{inv}}$ in Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021). The three components of $\boldsymbol {F}_{\text {NL}}$ are plotted in figure 3(d–f). The forcing $\boldsymbol {F}_{\text {NL}}$ gathers all nonlinear interactions of the forced flow with itself, including the Kelvin mode at order $A$, but also the particular shear flow solution at a lower order, $O(\varepsilon )$, given by (2.2), and the viscous correction to the Kelvin mode, also at a lower order, $O(A{ {Ek}}^{1/2})$, obtained by a series expansion of the wavenumber $k$ given in (2.9). Figure 3(d–f) shows an excellent agreement with the numerical results of Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) plotted in their figure 6(g–i). The largest component of $\boldsymbol {F}_{\text {NL}}$ is the axial component and it exhibits a negative and a positive lobe with a maximum value of $3.6 \times 10^{-3}$, while it is of $5 \times 10^{-3}$ in the simulations. The radial component has a weak dependence on $z$ and the azimuthal component is of order $Ek^{1/2}$ weaker than the other components, except for $\boldsymbol {F}_{BL}$ in the boundary layers.

The azimuthal component of the axisymmetric forcing, $F_{az.}=\boldsymbol {F}\cdot \hat {\boldsymbol {\theta }}$, plays a particular role here because it induces a geostrophic flow at an order ${ {Ek}}^{-1/2}$ higher than itself. The reason is that $F_{az.}$ directly forces a radial flow through the $\theta$-projection of the Navier–Stokes equation, $u_{az.} = \frac {1}{2} F_{az.}$, which itself forces an axial flow through the incompressibility condition

(5.2)\begin{equation} w_{az.}(r,z) ={-}\int \frac{1}{r} \frac{\partial(r u_{az.})}{\partial r} \text{d} z. \end{equation}

To cancel the $z$-odd part of this flow on the end walls, a geostrophic flow is necessary such that $(\tilde {w}_0)_{az.} = - w_{az.}$ in $z=-h/2$. This axial flow then yields a geostrophic flow

(5.3)\begin{equation} (v_0)_{az.}(r) ={{Ek}}^{{-}1/2} \int_{{-}h/2}^{0} F_{az.}\,\text{d} z. \end{equation}

The $\theta$-component of the axisymmetric forcing $F_{az.}$ therefore induces a geostrophic flow at order $F_{az.}{ {Ek}}^{-1/2}$, which is found to be proportional to the $z$-average of $F_{az.}$. Here, we note that this reasoning allows us to show that the nonlinear interactions between two Kelvin modes of different axial wavenumbers do not force a geostrophic flow. For instance, consider a non-resonant Kelvin mode of real amplitude $A_i=O(\varepsilon )$ and wavenumber $k_i = (2i - 1) {\rm \pi}/h$ interacting with the resonant Kelvin mode of amplitude $A$ and wavenumber $k = (2j-1) {\rm \pi}/h$. The nonlinear interactions will yield a non-zero azimuthal forcing $F_{az.}$ with terms proportional to $\varepsilon \, \mbox {Im}(A) \cos (k\pm k_i) z$. However, using (5.3), we see that this forcing does not give rise to any geostrophic flow because the sum or difference of axial wavenumbers, $(k\pm k_i)$, will always be an integer multiple of $2{\rm \pi} /h$ and will have a zero average along $z$.

The azimuthal components of the different forcing terms are plotted in figures 3(b,e) and 4(a–d). It shows that $\boldsymbol {F}_k\cdot \hat {\boldsymbol {\theta }}$ and $\boldsymbol {F}_{\text {NL}}\cdot \hat {\boldsymbol {\theta }}$ are both negative (except near $r=1$ for $\boldsymbol {F}_{\text {NL}}\cdot \hat {\boldsymbol {\theta }}$), even if the components $\boldsymbol {F}_{shear}\cdot \hat {\boldsymbol {\theta }}$ and $\boldsymbol {F}_{BL} \cdot \hat {\boldsymbol {\theta }}$ have positive average along $z$.

Figure 4.

Azimuthal forcing terms of the different components of the geostrophic flow: (a) $\boldsymbol {F}_{shear} \cdot \hat {\boldsymbol {\theta }}$; (b) $\boldsymbol {F}_{2k} \cdot \hat {\boldsymbol {\theta }}$; (c) $\boldsymbol {F}_{BL} \cdot \hat {\boldsymbol {\theta }}$; (d) $\boldsymbol {F}_{\text {NL}} \cdot \hat {\boldsymbol {\theta }}$, with $\boldsymbol {F}_{\text {NL}} = \boldsymbol {F}_{shear} + \boldsymbol {F}_{2k} + \boldsymbol {F}_{BL}$. Parameters are the same as in figure 3.

5.2. Zonal flow

The geostrophic flow $(\boldsymbol {u}_0)_{k}$ produced in response to $\boldsymbol {F}_k$ is plotted in figure 5(a). Its azimuthal velocity is always negative with a minimum value equal to $-2.2 \times 10^{-3}$ in good agreement with the minimum value of approximately $-2.5 \times 10^{-3}$ found numerically by Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) (see their figure 6o). However, in the numerical simulations, the azimuthal velocity vanishes at $r=1$ in order to satisfy the no-slip boundary conditions whereas it remains negative in the present theory.

Figure 5.

Mean azimuthal flow generated by different sources. (a) Velocity $(\boldsymbol {u}_0)_{k}$ given by (4.3) and induced by the forcing $\boldsymbol {F}_k$. (b) Velocity $(\boldsymbol {u}_0)_{shear} + \boldsymbol {u}_{2k} + (\boldsymbol {u}_0)_{2k} + (\boldsymbol {u}_0)_{Ek} + (\boldsymbol {u}_0)_{BL} = \boldsymbol {u}_{2k} + \boldsymbol {u}_0 - (\boldsymbol {u}_0)_{k}$ given by (4.8), (4.13), (4.14) and (4.19) and induced by the nonlinear forcing term $\boldsymbol {F}_{\text {NL}}$. (c) Total azimuthal velocity $\boldsymbol {u}_{2k} + \boldsymbol {u}_0$ induced by all sources. Parameters are the same as in figures 3 and 4.

The azimuthal flow $(\boldsymbol {u}_0)_{shear} + \boldsymbol {u}_{2k} + (\boldsymbol {u}_0)_{2k} + (\boldsymbol {u}_0)_{Ek} + (\boldsymbol {u}_0)_{BL} = \boldsymbol {u}_{2k} + \boldsymbol {u}_0 - (\boldsymbol {u}_0)_{k}$ in response to the nonlinear forcing term $\boldsymbol {F}_{\text {NL}}$ is plotted in figure 5(b). It exhibits two negative lobes of azimuthal velocity with a minimal value equal to $-1.6 \times 10^{-3}$. This value can be compared with the value of $-2.5\times 10^{-3}$ found numerically by Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) (their figure 6s), and the difference can be explained by the 30 % difference in the Kelvin mode kinetic energy.

Adding all responses yields the zonal flow $\boldsymbol {u}_{2k} + \boldsymbol {u}_0$ (figure 5c). It is almost everywhere negative except near the lateral wall around the equator. The velocity distribution is quantitatively similar to that found numerically by Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) in their figure 6(k). We find a minimum value of $-3.3 \times 10^{-3}$, close to the numerical value of $-4.5 \times 10^{-3}$ with again a difference that can mainly be imputed to the 30 % difference in kinetic energy.

Averaging this total flow over the height of the cylinder filters out the non-geostrophic response $\boldsymbol {u}_{2k}$ leading to the total geostrophic flow $\boldsymbol {u}_0$. This geostrophic flow, plotted in figure 6 as a black solid line, is in excellent agreement with the numerical result of Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) plotted as a black dotted line. This geostrophic flow contains several terms, the Coriolis term $(\boldsymbol {u}_0)_{k}$ being exactly opposite to the shear term $(\boldsymbol {u}_0)_{shear}$. Apart from these two terms, the dominant contribution comes from the nonlinear coupling between the inviscid Kelvin mode and its viscous correction in the bulk $(\boldsymbol {u}_0)_{Ek}$, which is retrograde. This flow is weakly compensated by the geostrophic adaptation $(\boldsymbol {u}_0)_{2k}$ to the zonal flow $vit_{2k}$ and by the nonlinear coupling in the boundary layers $(\boldsymbol {u}_0)_{BL}$. We note here that $(\boldsymbol {u}_0)_{2k}$ and $(\boldsymbol {u}_0)_{BL}$ are approximately equal for the parameters chosen; this is merely incidental.

Figure 6.

Geostrophic flow generated by different sources. The components $(\boldsymbol {u}_0)_{k}$, $(\boldsymbol {u}_0)_{shear}$, $(\boldsymbol {u}_0)_{2k}$, $(\boldsymbol {u}_0)_{Ek}$ and $(\boldsymbol {u}_0)_{BL}$, given by (4.3), (4.8), (4.12), (4.15) and (4.19), are represented by solid lines. They sum up to form the total geostrophic flow $\boldsymbol {u}_0$ (black solid line). The velocity is normalised by the kinetic energy of the Kelvin mode $E$. The dotted black line show the numerical results of Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) for the same parameters (see their figure 10). Parameters are the same as in figures 3 and 5.

5.3. Dependence of the zonal flow on the aspect ratio $h$

We have shown that, in the steady case, the geostrophic flow is proportional to $|A|^{2}$. This was also found numerically by Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021): in their figure 11(c), they showed that the ratio between the kinetic energy of the azimuthal flow $E_0$ and the square of the kinetic energy of the Kelvin mode $E^{2}$ is nearly independent of the Ekman number and the tilt angle. Here, we show that $E_0$ is proportional to $|A|^{4}$, while $E$ given by (4.25) is proportional to $|A|^{2}$, which is consistent with the observation of Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021).

Although the ratio $E_0 / E^{2}$ is independent of the forcing amplitude $\varepsilon$ and Ekman number ${ {Ek}}$, it depends on the aspect ratio $h$. This is because the aspect ratio influences the resonant frequency $\omega _{n,j}$ and the radial wavenumber $\delta _{n,j}$. In figure 7, we compare the quantity $E_0 / E^{2}$ extracted from Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) with our prediction, when $h$ is varied.

Figure 7.

Ratio between the kinetic energy $E_0$ of the zonal flow $\boldsymbol {u}_{2k}+\boldsymbol {u}_0$ and the squared kinetic energy $E$ of the forced Kelvin mode. This ratio is plotted as a function of the aspect ratio $h$ for the $(1,1)$ resonance. The symbols correspond to the numerical results obtained by Albrecht et al. (Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021) for ${ {Ek}}=2\times 10^{-4}$. The solid line corresponds to the present theory.

The theoretical prediction is found to be in excellent agreement with these numerical results although the theory tends to slightly underestimate the axisymmetric flow in the limit of small and large aspect ratios. This may be due to the presence of the boundary layers on the cylinder sidewalls.