L Oruba

A M Soward

E Dormy

Spin-down in a rapidly rotating cylinder container with mixed rigid and stress-free boundary conditions

A comprehensive study of the classical linear spin-down of a constant density viscous fluid (kinematic viscosity ν) rotating rapidly (angular velocity Ω) inside an axisymmetric cylindrical container (radius L, height H) with rigid boundaries, that follows the instantaneous small change in the boundary angular velocity at small Ekman number E = ν/H 2 Ω 1, was provided by Greenspan & Howard (J. Fluid Mech., vol. 17, 1963, pp. 385-404). E 1/2 -Ekman layers form quickly triggering inertial waves together with the dominant spin-down of the quasi-geostrophic (QG) interior flow on the O(E -1/2 Ω -1 ) time-scale. On the longer lateral viscous diffusion time-scale O(L 2 /ν), the QG-flow responds to the E 1/3 -side-wall shear-layers. In our variant the side-wall and top boundaries are stress-free; a setup motivated by the study of isolated atmospheric structures, such as tropical cyclones, or tornadoes. Relative to the unbounded plane layer case, spin-down is reduced (enhanced) by the presence of a slippery (rigid) side-wall. This is evinced by the QG-angular velocity, ω , evolution on the O(L 2 /ν) time-scale: Spatially, ω increases (decreases) outwards from the axis for a slippery (rigid) side-wall; temporally, the longtime ( L 2 /ν) behaviour is dominated by an eigensolution with a decay rate slightly slower (faster) than that for an unbounded layer. In our slippery side-wall case, the E 1/2 × E 1/2 corner region that forms at the side-wall intersection with the rigid base is responsible for a ln E singularity within the E 1/3 -layer causing our asymptotics to apply only at values of E far smaller than can be reached by our Direct Numerical Simulation (DNS) of the entire spin-down process. Instead, we solve the E 1/3 -boundarylayer equations for given E numerically. Our hybrid asymptotic-numerical approach yields results in excellent agreement with our DNS.

Introduction

Intense nearly axisymmetric vortices often develop in geophysical flows: tornadoes, or hurricanes in the atmosphere, while in the ocean, the Sea Surface Height variability appears dominated by westward-propagating mesoscale eddies throughout most of the World Ocean [START_REF] Chelton | Global observations of nonlinear mesoscale eddies[END_REF]. In order to understand their characteristics, it is instructive to model such objects as isolated structures (see Persing et al. 2015, and references therein). Their natural embedding inside a cylindrical domain introduces the need for an artificial outer circular side-wall boundary. Though no-slip boundaries could 2 L. Oruba, A. M. Soward and E. Dormy be contemplated, they introduce additional friction and so nowadays other boundary conditions are usually considered. Whatever the model, boundary layers are important (see, e.g., [START_REF] Smith | Hurricane boundary-layer theory[END_REF]. [START_REF] Williams | Thermal convection in a rotating fluid annulus: Part 3. Suppression of the frictional constraint on lateral boundaries[END_REF] was the first to suggest removing some frictional constraints associated with the side-walls in the numerical models, in order to more faithfully mimic the unbounded physical domain. Viscous friction was then only retained at the bottom boundary, as in atmospheric flows. Since then, stress-free side-walls, at which the angular velocity gradient vanishes, have been used in several numerical studies (see, e.g., Read 1986a,b); while matching to an external azimuthal flow corresponding to an axial line vortex (no viscous force), is sometimes preferred (see, e.g., [START_REF] Montgomery | Axisymmetric spindown dynamics of hurricane-like vortices[END_REF].

Here, we investigate the effects of a stress-free side-wall on such flows. Specifically, we consider the classical linear spin-down problem (see, e.g., [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] and the review of [START_REF] Benton | Spin-up[END_REF]; for non-linear studies see, e.g., [START_REF] Wedemeyer | The unsteady flow within a spinning cylinder[END_REF]; [START_REF] Hyun | Numerical solutions for spin-up from rest in a cylinder[END_REF] and the review of [START_REF] Duck | Spin-up of homogeneous and stratified fluids[END_REF]) in a domain with no-slip boundary at the bottom but modified by the presence of stress-free side-wall and top boundaries. Like [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF], we find that during the early stages, the mainstream spin-down exterior to all boundary layers is characterised by an angular velocity that is spatially constant but decays exponentially in time just as in an infinite (unbounded) plane layer (see §1.1). However on a somewhat longer time-scale the quasi-geostrophic angular velocity develops a radial (non-constant) structure dependent on the jump condition across the Stewartson (1957) E 1/3 -side-wall shear-layer. In the rigid boundary case studied by [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF], that shear-layer is passive and the quasi-geostrophic angular velocity simply vanishes at the mainstream boundary. This may be interpreted as a retrograde torque at the outer boundary that enhances the decay of the angular velocity with increasing radius. In our stress-free boundary case, the torque at the true boundary certainly vanishes. Nevertheless, the return meridional flow in the E 1/3 -layer is necessarily opposite to that needed to achieve spin-down in the mainstream. So from that point of view it is perhaps not surprising to find that the outer boundary essentially opposes spin-down rather than enhancing it. Essentially the E 1/3 -layer is no longer passive but completely alters the effective mainstream boundary condition. Consequently, the mainstream flow experiences a prograde boundary-torque proportional to the quasi-geostrophic angular velocity there. In turn, that causes the angular velocity to increase (rather than decrease) outwards (see (1.8)). This was an unexpected finding; though with hindsight it is easily explained.

The Greenspan and Howard problem

Here we describe in more detail our variant of the [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] model. Relative to cylindrical polar coordinates, (r , θ , z ), we consider a cylindrical container of height H and radius L rotating rapidly with angular velocity Ω about its axis of symmetry. The container is filled with constant density fluid of viscosity ν, which initially at time t = 0 rotates rigidly with angular velocity (1 + ε)Ω, where 0 < ε 1. The top boundary (r < L, z = H) and the side-wall (r = L, 0 < z < H) are impermeable and stress-free. The lower boundary (r < L, z = 0) is rigid. For that reason alone the initial state of rigid rotation (1 + ε)Ω of the fluid cannot persist and the fluid spins down to the final state of rigid rotation Ω of the container as t → ∞.

The rapid rotation of the system is measured by the small Ekman number

E = ν H 2 Ω 1 .
(1.1)

A spin-down problem 3 On very short time scales shear-layers form adjacent to the all boundaries of width δ (t ) = (νt ) 1/2 = H(EΩt ) 1/2 (1.2) due to viscous diffusion. On the short inertial wave time-scale t * = 1/Ω, a quasi-steady Ekman layer of width δ E = (ν/Ω) 1/2 = E 1/2 H forms adjacent to the lower boundary z = 0. Various inertial waves and quasi-geostrophic (z -independent; QG) motions exterior to all boundary layers are generated in the mainstream which decay primarily due to the ensuing Ekman suction into (or blowing out of) the Ekman layer. Our numerical results generally show that motion is soon dominated by the QG-flow which spins down on the longer time-scale H (νΩ) 1/2 = E -1/2 /Ω. Indeed, relative to the frame rotating with angular velocity Ω, the cylindrical components of the mainstream velocity are simply ȗ , v , w = εΩκ 1 2 σE 1/2 r , r , σE 1/2 H -z exp -σE 1/2 Ωt , (1.3a)

as determined by the solution to the initial value problem for a layer of unbounded horizontal extent outlined in appendix A and previously considered by [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF]. The purpose of appendix A is to extend Greenspan & Howards' κ ≈ 1, σ ≈ 1 results to determine their more precise forms (A 5b) and (A 8a) respectively. They have expansions

κ = 1 + 1 4 E 1/2 + O(E) , σ = 1 + 3 4 E 1/2 + O(E) (1.3b,c)
(see (A 9a) and (A 8b) respectively), which prove useful when comparing our numerical results at finite E with our asymptotic predictions. The spin-down motion (1.3a) is characterised by a fluid flux inside the Ekman layer towards the axis which is blown out into the mainstream with velocity w z =0 = ΩκσE 1/2 H. This generates a radial outflow in the mainstream of magnitude ȗ z >0 = 1 2 ΩκσE 1/2 Hr , which by conservation of angular momentum causes the angular velocity decay ω (t ) = v /r = εΩκ exp -σE 1/2 Ωt .

(1.4)

The mainstream outflow is blocked by the outer boundary r = L, where a quasisteady [START_REF] Stewartson | On almost rigid rotations[END_REF] side-wall layer of width δ S = E 1/3 H forms (see [START_REF] Barcilon | Stewartson layers in transient rotating fluid flows[END_REF]) on the time-scale t * = E -1/3 /Ω, long compared to the Ekman layer time-scale 1/Ω but short compared to the spin-down time E -1/2 /Ω. Following the formation of the E 1/3 -Stewartson layer, an ever thickening QG-shear-layer emerges of width δ (t ) = H(EΩt ) 1/2 (see (1.2)). On the spin-down time-scale E -1/2 /Ω, it has width E 1/4 H, which prompted [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] to refer to it as an E 1/4 -layer. The implied link with the static E 1/4 -Stewartson layer is misleading as the dynamical balances are different and there are no persistent (quasi-static) E 1/4 -layers in the spin-down problem (see also the discussions below (2.13) and in §5).

Following the establishment of the E 1/3 -layer (i.e., E 1/3 Ωt 1), the QG-shear-layer thickens: δ (t ) = H(EΩt ) 1/2 . This modifies the spin-down profile (1.4) both spatially and temporally on the longer lateral diffusion time

L 2 /ν = 2 E -1 /Ω, where = L/H (1.5)
is the container aspect ratio. This suggests a two time-scale approach (see, e.g., Hinch 1991, Chapter 7) in which (1.4) is modulated by a factor ω r /H, EΩt so taking the form ω (r , t ) = εΩκ ω r /H, EΩt exp -σE 1/2 Ωt .

(1.6)

Matching with the short time solution (1.4), as EΩt ↓ 0, is achieved by demanding that L. Oruba, A. M. Soward and E. Dormy initially ω(r /H, 0) = 1 .

(1.7)

The outer r = L boundary condition on the QG-mainstream flow depends on the jump conditions across the E 1/3 -side-wall layer.

For the rigid outer boundary considered by [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] the jump in azimuthal velocity across the E 1/3 -layer is negligible and so the outer boundary condition on the mainstream angular velocity is simply ω L, t = 0. It causes the early rigid rotation profile (1.4) to be further eroded to zero through lateral diffusion inwards from the outer boundary. For our stress-free outer boundary, consideration of the E 1/3 -layer shows that the boundary condition on the mainstream angular velocity is different of mixed type taking the Robin boundary condition form

∂ ω ∂r = α ω H at r = L , (1.8)
where α depends on the nature of the E 1/3 -layer and is a function of E alone, and 2πHL 3 (ρν)∂ ω ∂r r =L (viscosity ρν) is the total outer (mainstream) boundary torque on the QG-mainstream flow (see also (3.4)). The long time behaviour of the solution has the functional form

ω ≈ ω r /H exp k 2 -2 EΩt , (1.9) as -2 EΩt → ∞. Unlike α in (1.8), the constant k is a function of both E (via α) and .
The main point that we wish to stress is that whereas (1.3) defines the mainstream eigensolution on the spin-down time E -1/2 /Ω with growth rate -σE 1/2 Ω, that solution evolves over the lateral diffusion time-scale 2 E -1 /Ω and asymptotes to the form

ω ≈ εΩκ ω r /H exp -qEΩt , q = σE -1/2 -k 2 -2 , (1.10a,b)
with the different growth rate -qEΩ and non-constant radial profile proportional to ω r /H . Remarkably the stress-free outer boundary actually hinders the spin-down as evinced by the positive growth rate k 2 -2 EΩ of ω (see (1.9)). The phenomenon can be traced to the mainstream boundary torque proportional to ∂ ω ∂r r =L determined by (1.8) with α > 0. This long time behaviour is very different to that for the case of a rigid outer boundary for which the corresponding form for ω r /H, EΩt decays (see the n = 1 term of (3.8a)) rather than grow exponentially. The striking contrast made in our discussion §5 between (5.1) and (5.2), which incorporate these results, highlights the central theme of our paper.

Outline

Our notation is complicated and best outlined in general terms at the outset. Our basic variables, like the angular velocity are identified by the breve accent ˘. Only variables with a superscript star are dimensional; otherwise they are dimensionless. Azimuthal QG-flows independent of z * are identified by an overline . All these notations were introduced in §1.1. There the mathring accent ˚was also introduced in the restricted sense of mainstream QG-flows, being the amplitude that remains after the spin-down decay factor exp -σE 1/2 Ωt is removed (see (1.6)). Below, we will however use it in the general sense (see (2.7a)), which includes the mainstream QG-part together with the Ekman layer contribution denoted by the tilde (see (2.10)). The final eigensolution that emerges as t → ∞, like the QG-form (1.10) with the decay factor exp -qEΩt , is otherwise without accents (see (2.29)).

We formulate the transient flow problem in its entirety in §2. We outline the basic properties of the Ekman layer together with its consequences for the QG-flow in §2.1, particularly the QG-momentum equation (2.13) governing ω and its relation (2.20a) to the z -average ω . They differ by a small amount due to the angular momentum deficit in the Ekman layer. Though the difference is small, we find that our comparisons with the numerics and asymptotics is improved by use of (2.20a). We consider the eigenvalue problem for the final decay of the QG-flow (1.10a) in §2.2 from an asymptotic point of view on the basis that the eigenvalue k in (1.10b), related to α in (1.8), is known. Then in §2.3 we report results from direct numerical simulations, which confirm that the solution approaches the proposed asymptotic form. In §2.4, we identify the shape and decay rate but not the amplitude of the eigensolution again using numerical methods.

The transient evolution from the classical spin-down mainstream flow (1.3) to the final decay mode (1.10) (or more specifically (5.1)) is of particular interest. So in §3 we solve the angular momentum equation (2.13) for ω subject to the initial condition (1.7) on the basis that α in the boundary condition (1.8) is known. The numerical solution is described in §3.1, while a power series solution valid for short times is described in §3.2 which notably is valid on the spin-down time E -1/2 /Ω when the shear-layer width is E 1/4 H. On that time it provides our analogue of the rigid boundary solution (Greenspan & Howard 1963, eq. (6.3)). Indeed our power series is informative on much of the lateral diffusion time scale 2 E -1 /Ω and even provides a good approximation to ω at the outer boundary until the final asymptotic behaviour is clearly evident.

For a complete understanding of the transient evolution together with the final decay mode, we need a theory for α. To that end, in §4 we formulate the boundary layer equations for the E 1/3 -layer correct not simply to leading order but correct to O(E 1/6 ). We study the leading order problem asymptotically in §4.1.1 in the spirit of [START_REF] Stewartson | On almost rigid rotations[END_REF]. As the solution is singular near the bottom outside corner at (r , z ) = (L, 0), there it is best expressed in a similarity form which we provide and expand upon in §4.1.2. Unfortunately this leading order solution contains no QG-part. To extract the QG-part, we continue in §4.1.3 to O(E 1/6 ) at which the singularity forces a ln E type dependence in the asymptotic value (4.27) of α. As the numerical value of E needed for the applicability of the asymptotic theory is minute, we cannot compare directly asymptotic and numerical results. To overcome this obstacle we study in §4.2 the numerical solution of the side-wall boundary layer equations (4.3), containing the parameter ≡ E 1/6 explicitly at fixed E. Their numerical solution determines α, which we find agrees with that predicted by the full numerical eigensolution. In summary, we are unable to solve numerically the complete governing equations (2.4) at values of E small enough to reach the true asymptotic limiting behaviour described in §4.1 which is applicable in the limit E → 0. However, our numerical solutions of the asymptotically derived boundary layer equations (4.3) in §4.2 together with the transient results of §3 lead to a comprehensive and consistent picture of the asymptotic behaviour of the solution at small but finite E.

The mathematical problem

Whereas our primary concern is with the evolution of ω on the long lateral diffusion time-scale 2 E -1 /Ω together with the final eigenmode (1.10) state, we will nondimensionalise our governing equations on the short inertial wave time-scale 1/Ω as that is the time-scale on which inertial waves are manifest in the transient solution. Our numerical solution needs to take that evolution into account. Accordingly on writing (r , θ , z ) = (Hr, θ, Hz) , ȗ = εHΩ ȗ , t = t/Ω , (2.1a-c)

L. Oruba, A. M. Soward and E. Dormy the equations of motion are

∂ ȗ/∂t + 2 z × ȗ = ∇p + E∇ 2 ȗ , ∇ • ȗ = 0 , (2.2a,b)
where z is the unit vector in the z-direction (Ω ≡ Ω z) and p is a suitably nondimensionalised measure of pressure. The corresponding vorticity equation is

∂(∇ × ȗ)/∂t -2 z • ∇ ȗ = E∇ 2 (∇ × ȗ) . (2.2c) We write ȗ = ȗ , v , w = - E 1/2 r ∂ ψ ∂z , r ω , E 1/2 r ∂ ψ ∂r ,
(2.3a)

∇ × ȗ = -r ∂ ω ∂z , E 1/2 γ , 1 r ∂ ∂r (r 2 ω) , (2.3b) where γ = -∇ 2 - 1 r 2 ψ r = - 1 r D ψ , D = r ∂ ∂r 1 r ∂ ∂r + ∂ 2 ∂z 2 . (2.3c,d)
The scaling of the streamfunction ψ by the factor of E 1/2 anticipates our primary interest in QG-flows driven by Ekman suction. The azimuthal components of the momentum and vorticity equations (2.2) give

r 2 ∂ ω ∂t -2 E 1/2 ∂ ψ ∂z = E D r 2 ω , (2.4a)
∂ ∂t D ψ + 2r 2 E 1/2 ∂ ω ∂z = ED 2 ψ .
(2.4b)

The initial conditions are ψ = 0 , ω = 1 everywhere at t = 0 .

(2.5)

For t > 0 the boundary conditions are

ψ = ∂ ω ∂r = ∂ w ∂r = 0 at r = 0 and (0 < z < 1) , (2.6a) ψ = ∂ ψ ∂z = ω = 0 at z = 0 (0 < r < ) , (2.6b) ψ = ∂ 2 ψ ∂z 2 = ∂ ω ∂z = 0 at z = 1 (0 < r < ) .
(2.6c)

In order to discuss the evolution on the longer lateral diffusion time 2 E -1 /Ω, we make the further change of variables

ȗ , ψ , ω = κ ů , ψ , ω (r, z, τ ) exp(-σE 1/2 t) , t = 2 E -1 τ , (2.7a,b)
where we have incorporated the amplitude change κ and decay rate σE 1/2 (see (1.3b,c) respectively), predicted by the spin-down solution described in appendix A to occur on the shorter time E -1/2 /Ω (see also (1.4)). The new variables ů, ψ and ω satisfy (2.2)-(2.6) as before but with

∂ ∂t → -2 E ∂ ∂τ -σE 1/2 .
(2.8) So, for example, (2.4a,b) become

-2 E 1/2 r 2 ∂ω ∂τ -σr 2 ω -2 ∂ ψ ∂z = E 1/2 D r 2 ω , (2.9a) -2 E 1/2 ∂ ∂τ D ψ -σD ψ + 2r 2 E ∂ω ∂z = E 1/2 D 2 ψ . (2.9b)
We decompose the velocity ů into its mainstream ů and Ekman layer ů parts:

ů = ů + ů , (2.10)
with similar decompositions for the other variables. The QG-nature of the mainstream flow is dictated by (2.9b), which with (2.9a) determines

ω = ω(r, τ ) + O(E) , ψ = ψ(r, τ ) (z -1) + O(E) , ů = ů(r, τ ) + O(E 3/2 ) , (2.11a,b)
where, from (2.3a), (2.11b) and (2.9a) again,

-ψ = E -1/2 rů = 1 2 σr 2 ω -E 1/2 r 2 2 ∂ω ∂τ -D(r 2 ω) (2.11c)
and ω(r, τ ) is a function of r and τ yet to be determined. A fortunate feature of our two time-scale ansatz is that (2.11c) and (2.9b) are solved, correct to O(E 1/2 ) as advocated in (2.11a,b), by

-ψ = E -1/2 rů = 1 2 σr 2 ω , (2.12) ∂ω ∂τ = 2 r 2 D(r 2 ω) .
(2.13) Equation (2.13) highlights the absence of any E 1/4 -Stewartson layer, usually triggered by the Ekman suction term -ψ on the left-hand side of (2.11c). For our spin-down flow, that effect is exactly balanced by the inertial decay 1 2 σr 2 ω on the right-hand side (see (2.12)).

Ekman Layer

In the Ekman layer near z = 0 we set

v = rω(r, τ ) + v , ψ = ψ(r, τ ) (z -1) + ψ , ů = ů(r, τ ) + ů .
(2.14a,b)

Correct to leading order, the governing equation (2.2), as modified in (2.9), determines the Ekman layer equations

-E 1/2 σ v + 2 ů = E ∂ 2 v ∂z 2 , -E 1/2 σ ů -2 v = E ∂ 2 ů ∂z 2 ,                 r ů = -E 1/2 ∂ ψ ∂z , E 1/2 ψ + rů = -r z 0 ů dz , (2.15a,b)
where note has been taken of (2.11b), as well as E 1/2 ψ = E 1/2 ψ = -rů (see (2.12)) at z = 0. Equations (2.15a,b) must be solved subject to v = -rω, ů = -ů at z = 0 with v, L. Oruba, A. M. Soward and E. Dormy ψ, ů all tending to zero as z/E 1/2 ↑ ∞. We rewrite (2.15a) compactly as

E ∂ 2 Z ∂z 2 -2i -E 1/2 σ Z = 0 , Z = v -i ů , (2.16) with solution Z = -rω -iů exp -E -1/2 1 + 1 2 iE 1/2 σ 1/2 (1 + i)z .
(2.17)

From (2.15b) evaluated as E -1/2 z ↑ ∞ we obtain rů = -rRe ∞ 0 iZ dz = 1 2 E 1/2 r Re rω -iů (1 + i) 1 + 1 2 iE 1/2 σ 1/2 . (2.18)
Using (2.12), namely ů = 1 2 σE 1/2 rω, we recover the formula (A 8) for σ. Also of interest to us is the azimuthal Ekman layer flux

∞ 0 r ω dz = Re ∞ 0 Z dz = -1 2 E 1/2 rω + O(E) . (2.19) Using (2.14a) it means that the z-average of ω(r, z, τ ) is ω = µ + O(E) ω(r, τ ) , µ = 1 -1 2 E 1/2 + O(E) (2.20a,b) (see (A 9b))
, where

• • • ≡ 1 0 • • • dz (2.20c)
and µ has the specific definition (A 6b). So, in the case of ω = const. considered in appendix A, the relation ω = µω holds exactly without the O(E) correction mentioned in (2.20a). We stress this subtle difference as our Ekman layer calculation does not consider the O(E) Ekman layer corrections due to the effect of the term E(∇ 2 -∂ 2 /∂z 2 ) ů ignored in (2.15a,b). Whereas ω is the natural quantity to consider from the point of view of asymptotics, only ω can be measured unambiguously from our direct numerical simulation (DNS) of the complete governing equations (2.4) at finite E. For the largest value E = 3 × 10 -4 used, the small O E 1/2 difference is not insignificant and taking the correction into account improves the accuracy of our comparison of the DNS results with the asymptotics.

The final decay of the QG-flow: τ → ∞

As pointed out earlier, the r = mainstream boundary condition (1.8) only emerges after proper consideration of the E 1/3 -side-wall boundary layer, undertaken in §4 below.

Here we simply note that the final decay mode is described by (1.10a) which relative to the spin-down decay grows and takes the form (1.9):

ψ , ω (r, τ ) ≈ ψ , ω (r) exp k 2 τ as τ → ∞ . (2.21) It satisfies (2.13) when r d dr r d dr (rω) - kr 2 + 1 rω = 0 . (2.22)
In terms of the Modified Bessel function I 1 , the solution regular at r = 0 is

ω(r)/A 0 = ω0 (r) ≡ 2(kr/ ) -1 I 1 (kr/ ) , (2.23)
where the normalisation constant A 0 is chosen (see (2.25) and (2.27) below) such that ω0 (0) = 1. From (2.23) we deduce that (2.24c,d) where the prime denotes differentiation. Consideration of the transient QG-evolution in §3 below shows that

dω 0 dr ( ) = 2k I 1 (ρ) ρ ρ=k = 2I 2 (k) ≡ α ω0 ( ) , α = k I 2 (k) I 1 (k) , (2.24a,b) equivalently dω dr ( ) = α ω( ) , α = α/ ,
A 0 = 1 2 k 2 I 2 (k) k 0 ρ I 1 (ρ) 2 dρ k 2 I 2 (k) = k 0 ρ 2 I 1 (ρ) dρ (2.25) (see (3.7)). Integration of the identity 2ρ I 1 (ρ) 2 = -ρI 1 (ρ) 2 + 1 + ρ 2 I 1 (ρ) 2 (2.26a)
and noting again that ρI

1 (ρ) = ρI 2 (ρ) + I 1 (ρ) determines 2 k 0 ρ I 1 (ρ) 2 dρ = k 2 -2α -α2 I 1 (k) 2 (2.26b)
on use of (2.24b). Substitution of (2.26b) into (2.25) evaluates A 0 and whence from (2.23) we determine the boundary values

ω(r, τ ) exp -k 2 τ →        A 0 ≡ kα k 2 -2α -α2 1 I 1 (k) at r = 0 , A ≡ 2α k 2 -2α -α2 at r = .
(2.27)

Essentially in (2.23) we have identified the first term of a Dini series modal expansion (see [START_REF] Erdélyi | Higher Transcendetal Functions[END_REF] §7.10.4, especially eq. ( 49), and [START_REF] Watson | A Treatise on the Theory of Bessel Functions[END_REF] chapter XVIII). Finally to complete the solution (2.21), we note that (2.12) determines

-ψ = E -1/2 ru = 1 2 σr 2 ω = 1 2 1 + 3 4 E 1/2 + O(E) r 2 ω .
(2.28) 2.3. Direct numerical simulation (DNS) for the case = 1, E = 10 -4

We performed DNS of the governing equations (2.4) subject to the initial conditions (2.5) and boundary conditions (2.6). We will distinguish such solutions by the subscript DNS, i.e., our DNS-solution is ȗ , ψ , ω = ȗDNS , ψDNS , ωDNS (r, z, t). We solved (2.4) using second-order finite differences in space, and an implicit second-order backward differentiation (BDF2) in time. We used a stretched grid, staggered in the z-direction. Each simulation was initialised with a uniform distribution of ωDNS . The spatial resolution depends on the value of E, and was varied up to 2000 × 2000 to ensure convergence.

As τ → ∞ the final eigensolution takes the form outflow caused by Ekman blowing is revealed by the tilted ψ DNS -contours in figure 1(b), whose nature is consistent with the mainstream asymptotic result

κ -1 ȗDNS , ψDNS , ωDNS = u DNS , ψ DNS , ω DNS (r, z) exp(-q DNS 2 τ ) , ( 2
ψ = (z -1)ψ(r) ≈ -1 2 σ(z -1)r 2 ω(r) (2.30)
determined by (2.11b) and (2.28). This outflow is returned inside the E 1/3 -side-wall layer (see §4) , evident near r = ( = 1), to a sink, namely the E 1/2 × E 1/2 corner region in the neighbourhood of (r, z) = ( , 0). In turn that influx is ejected within the E 1/2 Ekman layer at the base of the E 1/3 -layer to provide the Ekman layer flux needed to adjust the QG-flow.

Numerical eigensolution (ENS) for the case = 1 for various E

As the DNS eigensolution (2.29) requires very considerable computational time, we found it more efficient to extract the functional form without detailed consideration of the transient evolution. By that expedient we were able to solve the eigenvalue problem for u (but only up to an arbitrary constant) and q (i.e., u DNS and q DNS in (2.29)) for the case = 1 for various values of E. We refer to it as the eigen (numerical) solution (ENS), for which [u , ψ , ω] = u ENS , ψ ENS , ω ENS (r, z) and q = q ENS . To place the results in context we note that our asymptotic theory predicts

ω (r, t) ≈ µω(r, t) ≈ (µκ) ω(r) exp(-q 2 τ ) as τ → ∞ (2.31) (see (2.7
), (2.20), (2.21) and (2.29)), where ω(r) is defined by (2.23). Since the notion of a mainstream QG-flow and an Ekman layer is an asymptotic (E 1) concept, ω ENS (r) evaluated at finite E is not clearly defined. To overcome this obstacle we simply set

ω ENS (r) ≡ µ -1 ω ENS (r)
(2.32) so that comparisons can be made with asymptotic theory for which q in (2.31) is unfortunately as yet unknown. However, on assuming that q = q ENS we may infer that the value of k characterising our analytic solution ω(r) (see (2.23)) is given by the upper boundary value ω ENS (r, 1) to equal ω(r). For that reason we compare the plots of ω(r), ω ENS (r) and ω ENS (r, 1) in figures 2(a-d) for various values of E. The plots confirm our mainstream expectations: ω(r) ≈ ω ENS (r) ≈ ω ENS (r, 1). Their values only differ significantly within the E 1/3 -side-wall (r = ) layer. The strength of that layer is indicated by the amplitude of the oscillations ω ENS (r, 1)-ω(r) = O E 1/6 . The magnitude of the oscillations ω ENS (r) -ω(r) = O E 1/3 are markedly smaller. The relative sizes are clearly visible on the plots, while the theory behind the orders of magnitude alluded to will be explained in §4.

k 2 = 2 σE -1/2 -q ENS (
We emphasise that so far the decay rate q has been predicted by the ENS. To obtain a closed form analytic solution we need to determine the value of α that characterises the mainstream boundary condition dω/dr( ) = αω( ) (see (2.24c)) by consideration of the E 1/3 -side-wall layer. With α known, the τ -growth rate k 2 (see (2.21)) is given by the solution k = k(α) of (k/ )I 2 (k)/I 1 (k) = α = α (see (2.24b,d)). That strategy, attempted in §4.1, is hindered by the asymptotic orderings ω(r, z) -ω(r) = O E 1/6 , µ -1 ω (r)ω(r) = O E 1/3 hinted at by figures 2(a-d). In short, though we are able to solve the vital O E 1/3 problem for the relatively small shear-layer correction µ -1 ω -ω, the tiny size of E needed for its validity is unreachable by our ENS. To bypass this difficulty, in §4.2 L. Oruba, A. M. Soward and E. Dormy we consider the shear-layer problem numerically retaining the relevant terms involving E needed to encompass both the O E 1/6 and O E 1/3 problems simultaneously. From this hybrid asymptotic-numerical method we determine α and in turn α = α and k = k(α) for various values of E. The ENS value q ENS and the hybrid asymptotic-numerical value q = σE -1/2 -k 2 -2 agree so well that for the purpose of plotting the graphs in figures 2(ad), they are essentially the same. The weak dependence of k and α on E evident from their values itemised in each respective figure proves to be a delicate issue that we discuss in §4, but particularly §4.1.

The transient evolution of the QG-azimuthal angular velocity

We consider the evolution of ω(r, τ ), which solves

∂ω ∂τ = 2 ∂ 2 ω ∂r 2 + 3 r ∂ω ∂r = 2 r 3 ∂Γ ∂r , Γ = r 3 ∂ω ∂r (3.1a,b) (see (2.13)), on the τ = O(1) time-scale.
Here ω(r, τ ) defines the actual QG angular velocity

ω(r, t) = κ ω(r, τ ) exp(-σE 1/2 t) (3.2) (see (2.7a))
, while Γ provides a measure of the total viscous couple on cylinders of fluid radius r. We need to solve (3.1) subject to the (τ = 0) initial condition

ω = 1 on 0 < r < , (3.3)
determined by the spin-down solution (1.3a), and the boundary conditions

1 2 Γ = ∂ω ∂r = 0 at r = 0 , α ω at r = . (3.4)
With α = α > 0 the outer boundary condition Γ ( , τ ) = 2 αω( , τ ) (see (2.24a,b)) says that azimuthal motion at the outer boundary leads to a couple that accelerates rather than brake the total rotational angular momentum r 3 ω on cylinders radius r. [START_REF] Park | Transient Stewartson layers of a rotating compressible fluid[END_REF] in their discussion of compressible Stewartson layers provided the complete transient solution of a diffusion equation similar to (3.1). For our purposes it is sufficient to note that, as τ → ∞, the asymptotic solution is given by (2.21). It is simply the first term of a modal expansion

ω(r, τ ) = A 0 ω0 (r) exp k 2 τ + ∞ n=1 A n ωn (r) exp -k 2 n τ , (3.5)
where ω0 (r) defined by (2.23) solves (2.22). Each ωn (r) (n > 0) satisfies the boundary conditions (3.4) and solves (2.22) with k 2 replaced by appropriate eigenvalues -k 2 n . It is easy to show that the corresponding eigenfunctions ωn (r) are orthogonal:

0 ωn ωm r 3 dr = 0 m = n . (3.6)
The initial condition ω(r, 0) = 1 thus determines

0 ωr 3 dr = A 0 0 ω 2 r 3 dr , (3.7)
which on substitution of ω0 (r) = 2(kr/ ) -1 I 1 (kr/ ) (see (2.23)) gives (2.25).

The solution corresponding to (3.5) for the [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] problem with the rigid outer wall, at which ω( , τ [START_REF] Park | Transient Stewartson layers of a rotating compressible fluid[END_REF], eq. ( 19)), where j n (> 0) are the zeros of the Bessel function J 1 .

) = 0 is ω(r, τ ) = ∞ n=1 2B n (j n r/ ) -1 J 1 (j n r/ ) exp -j 2 n τ (3.8a) (cf.
The slowest decaying mode n = 1 has

j 1 = 3.8317 • • • , (3.8b)
while its amplitude B 1 is determined in the same way as A 0 above and is given in (5.2).

3.1. Numerical results for the case = 1, E = 10 -4

We solved the reduced initial value problem (3.1), (3.3) and (3.4) for ω(r, τ ) numerically and compared the results with the z-average ωDNS (r, τ ) obtained from the DNS solution ωDNS (r, z, t) of the entire problem (2.4)-(2.6) for ω(r, z, t). To make that comparison we need to build on our analytic predictions. Firstly, the analysis of appendix A indicates that the amplitude of ω is modified by the factor κ (see (1.3a,b), also (A 9a)) on the spin-down time t = O E -1/2 , equivalently τ = O E 1/2 . Accordingly, we account for the spin-down decay rate and that amplitude modification in our analytic representation of the QG-mainstream angular velocity ω(r, t) in (3.2). Secondly, our asymptotic results predict ω = µ -1 ω (see (2.20a,b), also (A 9b)), and so we define the corresponding DNS value by ωDNS (r, t) ≡ µ -1 ωDNS (r, t)

(3.9) as in (2.32). Thirdly, in view of the definition ω = κ -1 ω exp σE 1/2 t (see (2.7)), we also define ωDNS (r, τ ) ≡ (µκ) -1 ωDNS (r, t) exp(σE 1/2 t)

for τ = O(1) , (3.10a)
where

(µκ) -1 = 1 + 1 4 E 1/2 + O(E) (3.10b) (see (A 9a,b)).
We illustrate the transient development of ω(r, τ ) and ωDNS (r, τ ) for aspect ratio = 1 and E = 10 -4 in figure 3 with other data itemised in the caption. Since (3.5) determines the late time behaviour

ω(r, τ ) ≈ A 0 ω0 (r) exp k 2 τ as τ → ∞ , (3.11)
we remove this exponential growth and plot ω(r, τ ) exp -k 2 τ and ωDNS (r, τ ) exp -k 2 τ instead at various times τ . The O(E 1/2 ) correction to ωDNS (see (3.10a)) that ensues through the factor (µκ) -1 (see (3.10b)) leads to very good agreement between the respective curves at each τ . Indeed, they are largely indistinguishable except in the E 1/3sidewall-layer, where ωDNS needs to adjust to meet the boundary condition ∂ω DNS ∂r = 0 at r = (= 1). Similar discrepencies are visible in figure 2(b). As noted in the early time asymptotics of the following §3.2 the left-hand (or centre) value ω(0, τ ) remains at unity until it eventually increases in response to diffusion across the domain driven by the right-hand (or outer) boundary condition. Eventually ω(0, τ ) grows exponentially (see (3.11)). To illustrate this behaviour, we plot log ω(0, τ ) versus τ together with its large τ asymptote (recall that ω0 (0) = 1) in figure 3(Inset). It suggests a transition between small and large τ behaviour over roughly the range 0.05 τ 0.1. This view is supported by figure 3 value unity at τ = 0. As time proceeds the flatness of the profiles near r = 0 for τ = 0.01 and 0.03 has largely disappeared by τ = 0.11. Subsequently the left (r = 0) and right (r = ) end point values decrease and increase (respectively) monotonically approaching the respective values A 0 and A (see (2.27)) as τ → ∞.

A short time (τ 1) series expansion

The initial condition ω = 1 (see (3.3)) satisfies the governing equation (3.1) and the boundary condition (3.4) at r = 0 but not at r = . The discontinuity at r = is readily accommodated by the similarity solution giving The similarity solution (3.13) is the first term of an asymptotic power series solution

∂ω ∂r = α erfc ζ , where ζ = 1 -(r/ ) √ 4τ , ( 3
1 (ζ) = ∞ ζ erfc ζ dζ = 1 √ π exp -ζ 2 -ζ erfc ζ . ( 3
ω = 1 + α ∞ n=1 (4τ ) n/2 n (ζ) valid on 0 < ζ < (4τ ) -1/2 (3.14)
for τ 1. We outline the solution built around (B 4)-(B 5) in appendix B. As a diagnostic for comparison with our direct numerical solution of (3.1)-(3.4) we note that the end point value (B 7) at r = determined by the terms (B 6) and (B 8)-(B 11), up to n = 4, is

ω( , τ ) ≈ 1 + 2α √ π τ 1/2 + α 3 2 + α τ + α √ π 5 2 + 4α + 4 3 α2 τ 3/2 + 3α 4 5 4 + 4α + 3α 2 + 2 3 α3 τ 2 + O(τ 5/2 ) . (3.15)
To compare with the large τ solution ω(r, τ ) = ω(r) exp(k 2 τ ), we expand ω(r, τ ) exp(-k 2 τ ) at both r = 0 and r = to the same order of accuracy:

ω(0, τ ) exp(-k 2 τ ) ≈ exp(-k 2 τ ) ≈ 1 -k 2 τ + 1 2 k 4 τ 2 + O(τ 3 ) , (3.16a) ω( , τ ) exp(-k 2 τ ) ≈ 1 + 2α √ π τ 1/2 + α 3 2 + α -k 2 τ + α √ π 5 2 + 4α + 4 3 α2 -2k 2 τ 3/2 + 3α 4 5 4 + 4α + 3α 2 + 2 3 α3 -k 2 α 3 2 + α + 1 2 k 4 τ 2 + O(τ 5/2 ) .
(3.16b)

Returning to results from the direct numerical solution of (3.1)-(3.4), whereas in figure 3(Inset) we plot the logarithm of the left-hand end-point value ω(0, τ ), in figure 4(a) we plot the right-hand end-point value ω( , τ ) exp -k 2 τ . We also show the approximations derived by using the power series representation of ω( , τ ) given by the series expansion (3.15) truncated at various levels. The first n = 1 term truncation (3.13) L. Oruba, A. M. Soward and E. Dormy identifies the τ 1/2 singular behaviour near τ = 0 but only gives a good approximation for very small τ . Further terms improve the τ -range of usefulness considerably. In figure 4(b), we employ the complete power series representation (3.16b) of ω( , τ ) exp -k 2 τ , which interestingly for even n (i.e., n = 2, 4) improves the approximation, exhibiting a longer τ -range of usefulness at each level of truncation. The improvement is not apparent for odd n (i.e., n = 1, 3). It is worth noting that the best (n = 4) truncation provides a reasonable approximation up until τ ∼ 0.1, by which time the solution portrayed in figure 3 is close to its final asymptotic τ → ∞ form. It must be appreciated that the power series solution (B 7) described in appendix B is truly asymptotic, not least because it takes no account of the left-hand (r = 0) boundary condition.

The E 1/3 -side-wall shear-layer

We pointed out in §1.1 that E 1/3 -shear-layer adjacent to the outer cylinder wall r = is established on the time scale t ∼ E -1/3 short compared to the spin-up time t ∼ E -1/2 . Once established its longer time evolution is best investigated via the modulation amplitudes ω(r, z, τ ) and ψ(r, z, τ ) of the spin-down solution, which satisfy the scaled governing equations (2.9a,b). Our ultimate objective is to establish analytically the value of α in (2.24c), that determines the mainstream boundary condition dω/dr = αω at r = . This relation is only meaningful for t E -1/3 by which time α is well defined and independent of t. So to determine α it is sufficient to consider the final time eigenfunction κ -1 ψ , ω exp q 2 τ = ψ , ω (r, z) . . Though formulated in this restricted sense for clarity, our shear layer analysis pertains to all t E -1/3 upon making the approximation q 2 τ = E 1/2 t in (4.1).

We restrict attention to the shear-layer flow outside the Ekman layer and represent it as the sum of the QG-mainstream solution (2.11a,b) and a shear-layer correction:

ω = ω(r) + ω( )V (ζ, z) , ψ = ψ(r) (z -1) + 2 ω( ) Ψ (ζ, z) , (4.2a,b)
where

ψ(r) = -1 2 σr 2 ω(r) , ζ = -2 (r -) and = E 1/6 (4.2c-e)
(see (2.23) and (2.28), also (2.30)). Here ζ is the stretched radial boundary layer coordinate and is the natural expansion parameter: E 1/3 = 2 and E 1/2 = 3 . The boundary layer corrections V , Ψ must tend to zero as ζ → -∞ as well ensure that the boundary conditions (2.6a) at r = are met.

To extract the boundary layer equations correct to O( ), it is sufficient to consider (2.9a,b) with σ ≈ 1 (so, e.g., (4.2c) becomes ψ(r) = -1 2 r 2 ω(r)). They determine

-2 ∂Ψ ∂z = ∂ 2 V ∂ζ 2 + V , 2 ∂V ∂z = ∂ 2 ∂ζ 2 ∂ 2 Ψ ∂ζ 2 + Ψ . (4.3a,b)
To the same order of accuracy, the boundary conditions (2.6a) are

Ψ = 1 2 (z -1) , ∂ 2 Ψ ∂ζ 2 = 0 , ∂V ∂ζ = -α at ζ = 0 , (4.4a) Ψ → 0 , V → 0 , as ζ → -∞ (4.4b
) for 0 < z < 1, where α = dω/dr( ) ω( ), while the Ekman jump condition at z = 0 on the flow exterior to the Ekman layer and the symmetry conditions (2.6c) at z = 1 are respectively

Ψ = 1 2 V at z = 0 , (4.4c) Ψ = ∂ 2 Ψ ∂z 2 = ∂V ∂z = 0 at z = 1 (4.4d ) for ζ < 0.
In view of the boundary conditions (4.4b), we may introduce the new variable

Θ(ζ, z) = ζ -∞ V dζ =⇒ V = ∂Θ ∂ζ , (4.5a,b)
and integrate (4.3b). By this device, we may symmeterise the shear-layer equations (4.3a,b) and express them in the compact complex form

2i ∂Ξ ∂z = ∂ ∂ζ ∂ 2 Ξ ∂ζ 2 + Ξ , (4.6a)
where

Ξ = Θ + i Ψ giving Υ = V + i W = ∂Ξ ∂ζ W = ∂Ψ ∂ζ (4.6b,c)
which is also useful.

Asymptotic approach

We show in §4.1.1 that, at zeroth order, the solution for V (ζ, z) has no z-mean, or more precisely, V = O( ). The objective of the first order problem investigated in §4.1.3 is to solve (4.25) for V which, in turn, determines our key quantity of interest α = --1 d V dζ ζ=0 , namely the z-average of (4.4a) 3 .

The zeroth order problem

On setting = 0, the solution of (4.3)-(4.6) is

Ξ(ζ, z) = - 2 √ 3π ∞ n=1 (-1) n n E k n ζ; -π/6 exp inπ(z -1) , (4.7a) Υ (ζ, z) = - 2 √ 3π ∞ n=1 (-1) n k n n E k n ζ; π/6 exp inπ(z -1) , (4.7b)
where

k n = (2nπ) 1/3 , E k n ζ; α = exp( 1 2 k n ζ) cos √ 3 2 k n ζ + α . (4.7c,d)
Verification of the essential relation Υ = ∂Ξ/∂ζ (see (4.6c)) follows from the property

d dζ E k n ζ; α = k n E k n ζ; α + π/3 =⇒ d 3 E dζ 3 = -k 3 n E (4.8a,b)
also. Since Im{exp[inπ(z -1)]} = 0 at z = 0 and 1, the boundary conditions (4.4c), namely Ψ (ζ, 0) = 0, at the bottom and (4.4d) at the top are all obviously met. Another property of the solution on z = 0, useful later in §4.1.

3 (see (4.26b)), is ζ -∞ Θ(ζ, 0) dζ = ζ -∞ Ξ(ζ, 0) dζ = - 2 √ 3π ∞ n=1 1 nk n E k n ζ; -π/2 , (4.9)
where we have used (4.8a). On use of (4.8b) it is readily seen that Ξ, defined by (4.7a), solves (4.6a). Furthermore, since d

2 dζ 2 E k n ζ; -π/6 ζ=0 = k 2 n E k n ζ; π/2 ζ=0 = 0 , (4.10)
it is clear that the two boundary conditions ∂ 2 Ψ/∂ζ 2 = ∂V /∂ζ = 0 at ζ = 0 are met, i.e.,

∂ 2 Ξ ∂ζ 2 (0, z) = ∂Υ ∂ζ (0, z) = 0 . (4.11) Furthermore, (4.7a,b) evaluated at ζ = 0 determines Ξ(0, z) = - 1 π ∞ n=1 (-1) n n exp inπ(z -1) = 1 π ln 1 + exp iπ(z -1) = π -1 ln 2 sin 1 2 πz + 1 2 i(z -1) (4.12a) = 1 π ln πz + i 2 (-1 + z) + ∞ n=1 Ξ n z 2n , (4.12b) ∂Ξ ∂ζ (0, z) = Υ (0, z) = - 2 1/3 π 2/3 ∞ n=1 exp(inπz) n 2/3 = - 2 1/3 π 2/3 Li 2/3 e iπz = - 2 1/3 Γ (1/3)e iπ/6 πz 1/3 + ∞ n=0 Υ n z n , (4.12c)
in which we have introduced the Polylogarithm Li 2/3 (η) = ∞ n=1 n -2/3 z n , and where

Ξ n = - ζ(2n) 2 2n nπ , Υ n = - 2 1/3 π 2/3 ζ (2/3) -n (iπ) n n! , (4.12d,e) in which ζ(m) = ∞ n=1 n -m
, for m > 1 and by analytic continuation elsewhere is the Riemann-Zeta function. The result (4.12a) verifies that the boundary condition Ψ (0, z) = 1 2 (z -1) is met, while in addition determining the result Θ(0, z) = π -1 ln 2 sin 1 2 πz . The series representations (4.12b,c) follow on use of (http://dlmf.nist.gov/25.8.E8), (http://dlmf.nist.gov/25.12.E12) respectively (here and below, we use this style of online reference to equations in Abramowitz & Stegun 2010) and are valid over the entire range 0 < z 1.

The power series solution

The power series forms (4.12b,c) of the solution at ζ = 0 identify clearly a singularity at (ζ, z) = (0, 0). To resolve it, we reconstruct the full solution (4.7a,b) by solving the shear-layer equation (4.6a) subject to the boundary conditions (4.11) and (4.12b,c) at ζ = 0. Those suggest the expansion

Ξ(ζ, z) = 1 π ln(πz) - i 2 + Y(Φ) + z i 2 - 1 6 Φ 3 + Ξ(Φ, z) , (4.13a) Υ (ζ, z) = z -1/3 Z(Φ) - 1 2 z 2/3 Φ 2 + Υ (Φ, z) , (4.13b)
where Φ is the similarity variable

Φ = ζ z 1/3 . (4.13c)
Each of the three terms in (4.13a) when grouped with the corresponding term in (4.13b) solves (4.6a). The first terms, the similarity solution, resolve the singularity. The second and the third terms, namely the sums

Ξ N (Φ, z) = N n=0 Υ n Y + n (Φ)z n+1/3 + [N/2] n=1 Ξ n Y 2n (Φ)z 2n , (4.14a) Υ N (Φ, z) = N n=0 Υ n Y n (Φ)z n + [N/2] n=1 Ξ n Y - 2n (Φ)z 2n-1/3 (4.14b) ([N/2] denotes integer part of N/2) as N → ∞,
Y(0) = 0 , Y n (0) = 1 , Y ± n (0) = 0 . (4.15a-c)
It is our belief that, like the exponentials which they approximate, the series (4.14) have an infinite radius of convergence (in ζ) and hence meet the boundary condition Ξ → 0 as ζ → -∞.

The required generating polynomial

Y n (Φ) = n k=0 (2i) k n! (3k)!(n -k)! Φ 3k (4.16a,b) with Y n (0) = 1 (see (4.15b)) solves L n Y n = 0 , where L n ≡ d 3 dΦ 3 + 2 3 i Φ d dΦ -3n . (4.16c,d) It determines Y + n (Φ) = - 0 Φ Y n (Φ) dΦ , and
Y - n (Φ) = dY n dΦ , (4.17a,b)
which solve

L n±1/3 Y ± n = 0 , (4.17c)
and have the property Y ± n (0) = 0 (see (4.15c)). The similarity form π -1 ln(πz) -1 2 i + Y(Φ) solves (4.6a) with = 0 when

L 0 Y = 2i π ⇐⇒ d 2 Z dΦ 2 + 2 3 iΦZ = 2i π , (4.18a,b)
where the limits of integration in

Y(Φ) = - 0 Φ Z(Φ) dΦ (4.18c)
are chosen such that Y(0) = 0 (see (4.15a)). The solution of (4.18b), which satisfies the remaining boundary conditions dZ/dΦ(0) = 0 and

Z → 0 as Φ → -∞ , (4.19a,b) is Z(Φ) = - √ 3 2 3 1/3 e iπ/3 Hi i 2 3 1/3 Φ + Hi 2 3 1/3 e -iπ/6 Φ , (4.20a)
where

Hi(z) = 1 π ∞ 0 exp -1 3 t 3 + zt dt , with Hi(0) = Γ (1/3) 3 2/3 π , (4.20b,c)
is the Scorer function (http://dlmf.nist.gov/9.12.E20), which satisfies Hi -z Hi = 1/π , (4.20d)

where the prime denotes derivative. Together (4.20a,c) determine (4.21) as required by the boundary condition (4.12c).

Z(0) = - 2 1/3 Γ (1/3) π e iπ/6 ,
It is worth reflecting at this point, that the governing equation (4.18b) is essentially that given by [START_REF] Dormy | Mathematical aspects of natural dynamos in The fluid mechanics of astrophysics and geophysics[END_REF] eq. (3.26b) in their interpretation of the point source similarity solution derived by [START_REF] Moore | The structure of free vertical shear layers in a rotating fluid and the motion produced by a slowly rising body[END_REF]. However, our lowest order combined boundary conditions ∂V /∂ζ = 0 and Ψ = -1 2 on ζ = 0 break the symmetry of the classical [START_REF] Stewartson | On almost rigid rotations[END_REF] and [START_REF] Moore | The structure of free vertical shear layers in a rotating fluid and the motion produced by a slowly rising body[END_REF] solutions. So whereas, the [START_REF] Moore | The structure of free vertical shear layers in a rotating fluid and the motion produced by a slowly rising body[END_REF] solution only involves a single Scorer function (see Dormy & Soward 2007, eq. (3.26a)), we need two (see (4.20a) so complicating matters. From a more general point of view, our solutions are valid on the semi-infinite interval -∞ < Φ 0 and would diverge if extended to 0 Φ < ∞. By contrast, [START_REF] Moore | The structure of free vertical shear layers in a rotating fluid and the motion produced by a slowly rising body[END_REF] considered similarity solutions valid on the infinite interval -∞ < Φ < ∞. These have proved useful in other E 1/3 shear-layer contexts (see, e.g., Marcotte et al. 2016, eqs. (4.1)-( 4.3)).

Substitution of (4.20a) into the integral (4.18c) determines

Y(Φ) = √ 3 e -iπ/6 J -i 2 3 1/3
Φ + e iπ/6 J -2 3

1/3 e -iπ/6 Φ , (4.22a) in which

J(z) = z 0 Hi(-z) dz . (4.22b)
For large negative Φ the asymptotic result (http://dlmf.nist.gov/9.12.E31) determines

Y(Φ) ≈ 1 π 3 ln(-Φ) + 2γ + ln 2 + 1 √ 3 + i 2 + 3i π(-Φ) 3 (-Φ 1) , (4.23)
where γ = 0.57721566

• • • is Euler's constant.
Values close to the bottom boundary z = 0 are important. On substitution (4.23) into (4.13a), the real and imaginary parts determine

Θ z↓0 ≈ 1 π 3 ln(-ζ) + 2γ + ln(2π) + 1 √ 3 , Ψ z↓0 ≈ 0 (4.24a,b)
(note that (-Φ) 3 z↓0 = 0, provided ζ = 0), from which we also deduce that

V z↓0 = ∂Θ ∂ζ z↓0 ≈ 3 πζ for 1 -ζ z 1/3 . (4.24c)
Also of interest is the integral (4.9), which evidently vanishes at ζ = 0. Its approach to zero is determined by (4.24a), which determines

ζ -∞ Θ z↓0 dζ ≈ 3 π ζ ln(-ζ) for 1 -ζ z 1/3 , (4.24d,e)
needed for our understanding of V (ζ) (see (4.26b)) in the next §4.1.3). Of course, the limit z ↓ 0 must be interpreted in the sense of our solution being valid outside the Ekman layer. As the Ekman layer width is O(E 1/2 ), we expect (4.24c) to be valid for -ζ E 1/6 (equivalently -r E 1/2 ). In appendix C, we evaluate Ψ and V (see (4.13a,b)) using the power series (4.14a,b) (see (C 1a,b)) and compare our results with those from the Fourier series (4.7a,b) in figure 5 below. Though the Fourier series needs a very large number of terms (10 5 ) to give an accurate representation near (ζ, z) = (0, 0), only a very few terms (N = 4) of the power series are needed to give a very accurate solution there. Indeed it is remarkable how good the power series solution is throughout the domain illustrated in figure 5. On increasing -ζ serious discrepancies on the entire range 0 < z 1 only begin to emerge at about ζ = -2.0, by which distance the values of both Ψ and V are small and relatively close to zero. We find this high degree of accuracy from such a low truncation quite striking. Indeed our belief is that the differences visible are due to the truncation and that the Fourier series and power series solutions coincide in the limit N → ∞.

The streamfunction Ψ in figure 5(a) shows uniform blowing from the right-hand boundary 0 < z < 1, where Ψ = 1 2 (z -1) (see (4.4a) 1 ) which, because Ψ = 0 on z = 0 (see (4.4b)), is all returned into the bottom right-hand corner (ζ, z) = (0, 0) so forming half an eddy. This is followed by a reverse flow eddy also spawned at (ζ, z) = (0, 0). Further eddies follow but they are essentially too small to be visible. The V -contours in figure 5(b) reflect the response of V to the Coriolis acceleration caused by the radial flow, proportional to -∂Ψ/∂z. The singularity in the corner identified by both V ≈ 3/(πζ) on z = 0 (see (4.24c)) and V ∝ z -1/3 on ζ = 0 (see (4.12c) and the red dashed curve, restricted to small z, in figure 8(b) below) gives the distinctive corner structure. The vertical nature of the V -contours close to the boundary at moderate -ζ (≈ 2.0) over a thin z-width is suggestive of an Ekman layer. That is an illusion, as in our asymptotics the Ekman layer has zero width. Rather it reflects how small z needs to be to achieve the large -Φ asymptotic behaviour. Though not visible in the figure, this feature continues up to ζ = 0, where it is truly a characteristic of the similarity solution V = z -1/3 Re{Z(Φ)}, which determines the V ≈ 3/(πζ) behaviour alluded to.

The first order (QG) problem

As the zeroth order solution V (ζ, z) = Re{Υ (ζ, z)} defined by (4.7b) has no mean part ( V = 0), we need to consider the next O( ) problem. On use of the boundary conditions (4.4c,d), the z-mean of (4.3a) determines

V (ζ, 0) = d 2 V dζ 2 + V , where V = V . (4.25a,b)
Since our shear-layer solution only concerns the flow outside the Ekman layer, we use the bar-notation for consistency with our ω notation for the mainstream QG-solution.

As V = O( ), the additional factor is included in (4.25b) for convenience. So ignoring the term V in (4.25) on the basis that it is O( ), integration with respect to ζ, noting the definition (4.5a), determines a. as ζ ↑ 0, the divergence of which means that α is not defined. This failure of our zeroth order asympotics is limited to the corner region

dV dζ (ζ) = Θ(ζ, 0) =⇒ V (ζ) =
-r = O(E 1/2 ), z = O(E 1/2 ), equivalently -ζ = O( ), z = O( 3 ). It identifies the ζ-length O(
), at which we must cut-off ln(-ζ) to obtain ln( ) + O(1), and whence

α = α FS ≈ -Θ(-, 0) ≈ - 3 π ln + O(1) = 1 2π ln E -1 + O(1) , (4.27)
where the subscript FS identifies the Fourier series origin (4.7) of this result. Recall that (4.27) is obtained from the similarity solution (4.22a), which is part of the power series solution (4.13) believed to be equivalent to (4.7). Henceforth, all results based on either (4.7) or (4.13) will be likewise distinguished.

Hybrid asymptotic-numerical solutions

Though the result (4.27) is asymptotically sound, it is impossible to determine the O(1) contribution to α FS without solving the problem for the E 1/2 × E 1/2 corner region, which we have not attempted. Indeed, as = E 1/6 is our expansion parameter, it is quite clear that no agreement with the full ENS is to be expected. Not only is ≈ 0.1468 for E = 10 -5 (as employed in figure 2(d)) not small, certainly ln itself is an O(1) number like the error neglected. To bypass this difficulty, we solved the problem posed by (4.3) and (4.4) numerically at small but finite . As a preliminary test of the approximation, we undertook the numerical solution with = 0. Despite the singularity at the corner (ζ, z) = (0, 0), the numerical results elsewhere agreed well with the Fourier series results portrayed in figure 5. Certainly to graph plotting accuracy the contours of constant Ψ and V are indistinguishable and so have not been re-plotted.

Henceforth our results obtained from our blend of asymptotic theory and numerical simulations (ANS) will be labelled by the subscript ANS. Though our implementation of ANS, adopts the boundary condition (4.4a) 3 :

∂V ANS ∂ζ (0, z) = -α ANS , the value α ANS = - dV ANS dζ (0) (4.28a,b)
is part of the answer. Indeed this consideration highlights the complication that the = 0 results, portrayed in figure 5, do not determine α ANS . Furthermore, despite appearances from our scalings, the value of α ANS is itself a function of = E 1/6 (asymptotically -ln , see (4.27)).

The purpose of the ANS-problem is to retain all the leading order elements of the asymptotic FS-problem, but also include the important O( ) terms. To this end we considered in detail the aforementioned case E = 10 -5 . Even though ≈ 0.1468 is only moderately small, the comparisons made in figures 5-7 are really encouraging; remember that the errors in the ANS-method are of order 2 ≈ 0.0215, i.e., 2%.

Since V FS (i.e., V ANS , when E = 0) displayed in figure 5(b) has no mean part, we plot the fluctuating part

V ANS (ζ, z) = V ANS (ζ, z) -V ANS (ζ) (4.29a)
in figure 6(b). Their constant-V ANS contours compare well with those for V FS in figure 5(b), but their numerical values differ by amounts that increase as -ζ approaches zero. The notion of a fluctuating part of Ψ ANS is less clear. So noting that V ANS measures the departure from quasi-geostrophy, we choose to plot the quantity V ANS (ζ, 0) 1 2 (z -1) is the QG-part stemming from the Ekman suction boundary condition ζ, 0). Furthermore as the non-zero value of Ψ ANS (0, 0) alters the amount of fluid flux that the shear-layer carries, we accommodate that effect by incorporating the corresponding Fourier solution combination -Ψ FS (ζ, z)+ 1 2 (z -1) scaled by the lower boundary velocity V ANS (ζ, 0). The upshot is that Ψ ANS defined by (4.29b) satisfies the boundary conditions Ψ ANS (ζ, 0) = 0 and Ψ ANS (0, z) = 1 2 (z -1), exactly like the = 0 Fourier series Ψ FS = Im{Ξ} (see (4.6b)). Note too that Ψ ANS → Ψ FS and V ANS → V FS as → 0. The comparison of Ψ ANS , V ANS for non-zero (E = 10 -5 ) in figures 6(a,b) with Ψ FS , V FS in figures 5(a,b) is topologically very good. However, quantitative agreement is moderate, which is to be expected as the expansion parameter ≈ 0.1468 is not really small. It should be emphasised that, whereas V ANS is a natural quantity, Ψ ANS is not, but rather has been constructed to provide an analogue to V ANS .

Ψ ANS (ζ, z) = Ψ ANS (ζ, z) + V ANS (ζ, 0) -Ψ FS (ζ, z) + 1 2 (z -1) , ( 4
Ψ ANS (ζ, 0) = 1 2 V ANS (
Finally in figure 7 we compare the ANS results Ψ ANS , V ANS with the corresponding ENS-quantities V ENS and Ψ ENS , defined by

ω( ) V ENS (ζ, z) = ω ENS (r, z) -ω(r) , (4.30a) 2 ψ( ) Ψ ENS (ζ, z)(r, z) = ψ ENS (r, z) -ψ(r)(z -1) (4.30b)
(see (4.2a,b) and recall that ψ(r) = -1 2 σr 2 ω(r), (4.2c)), for the case = 1 and again E = 10 -5 . To extract the best possible approximation for Ψ ENS from ψ ENS , we use σ = 1 + (3/4)E 1/2 (see (1.3c)). The agreement is good and consistent with the 2% estimated error from the neglect of terms order O( 2 ) neglected in the ANS-method. This success of the hybrid asymptotic-numerical method permits us to obtain results for small E and overcomes the obstacle faced by the pure asymptotic method which requires E to be essentially infinitesimally small.

An obvious test that addresses the question "How useful is the zeroth order asymptotic approximation that gives V FS ?" is to compare its lower z = 0 and end ζ = 0 boundary values with V ANS for a particular value of E. We do that in figures 5(a,b) for E = 10 -5 ( ≈ 0.1468) and = 1, as in figures 6 and 7. More precisely we plot 2Ψ ANS (ζ, 0) = V ANS (ζ, 0) in figure 8(a), where the -scaling of V is suggested by the corner boundary value 2Ψ ANS (0, 0) = -1 imposed by the end wall (ζ = 0) boundary condition (4.4a) 1 . Evidently this terminal value is approached but not reached by the ANS results and reflects the difficulties encountered by the numerics close to the corner singularity. The

V ANS (ζ) = -1 V ANS , the FS-solution V FS (ζ) = ζ -∞ Θ FS (ζ,
0) dζ determined by (4.9) and (4.26b) and the ENS-solution V ENS (ζ) which together with its derivative dV ENS /dζ are defined just as in (4.30a,b) by

E 1/3 ω( ) V ENS (ζ) = ω ENS (r) -ω(r) , (4.31a) ω( ) dV ENS dζ (ζ) = dω ENS dr (r) - dω dr (r) , (4.31b)
where ω ENS (r) is defined by (2.32). Here, ω ENS (r) and ω(r) in (4.31a), for our case E = 10 -5 , are plotted in figure 2 The lowest order Fourier solutions V FS and dV FS /dζ are independent of E and the reason for the choice of scaling in figures 9(a,b). Very good agreement is obtained between the ANS and ENS results just as in figure 8(b). By contrast the Fourier series prediction is acceptable for ζ -1.5 but is only fair for -1.5 ζ < 0. This highlights the fact that V is a blow up by a factor -1 of the small quantity V . This magnifies any errors and in part explains why it is difficult to obtain numerically accurate results from the asymptotics.

Our main objective is to determine α = dω/dr( )/ω( ) (see (2.24c)). We compare the

ANS and ENS values

α ANS = - dV ANS dζ (0) and α ENS = - dV ENS dζ (0) ( = 1) (4.32a,b)
for various values of E on figure 10. The remarkably good agreement of α ANS and α ENS in figure 10 appears to reflect the apparently fortuitous agreement of these derivatives at ζ = 0 evident on figure 9(b) for the particular case E = 10 -5 . For neighbouring values of ζ the agreement is not quite so good. Nevertheless, the generally favourable agreement of the ANS and ENS results illustrated by figures 9(a,b) and 10 vindicate our use of the hybrid asymptotical-numerical method-ANS. The corresponding true asymptotic value α FS given by (4.27) can only be extracted though the limiting procedure described thereabouts. In fact, our derivation of (4.27) assumes that α FS is determined, correct to the leading logarithmic order (but not the additional O(1) constant part), by that appeared to fit the ANS and ENS results. This differs by a factor 1/2 from the asymptotic value α FS = -(2π) -1 ln E predicted by (4.27). It suggests, of course, that the power law (4.33) is an illusion representing a snapshot over a very limited range of which is not small anyway! One therefore concludes that α FS is the true asymptotic limit but can only be reached by numerical values of that are essentially infinitesimals. Ultimately all our difficulties can be traced, as we remarked at the beginning of this subsection, to the singularity at (ζ, z) = (0, 0) that is only properly resolved by consideration of the E 1/2 × E 1/2 corner region in its vicinity.

- dV FS dζ (ζ) ≈ - 3 π ln(-ζ) , (4.32c

Discussion

As we explained at the outset, our investigation was motivated by the need to understand the role of artificial lateral boundaries in the modelling of essentially unbounded atmospheric flows. Our investigation of the spin-down problem suggests that the introduction of a stress-free outer boundary is far from adequate to obtain realism because motion is seriously affected by the presence of any impermeable boundary. The point here is that the classic spin-down mode (1.3a) exhibits a radial QG-outflow u ∝ r . So any boundary condition that prevents that natural outflow will have a profound effect as our analysis has shown. Indeed spin-down is a feature of vortex line shortening (equivalently fluid columns broadening). However the forced return flow in the E 1/3 -side-wall shearlayers has exactly the opposite effect causing a slight slow down of the spin-down process which is more pronounced towards the outer boundary. In summary good modelling of the unbounded system needs careful consideration of the fluid flux side-wall boundary condition upon which some results will be sensitive.

We have studied the entire transient process numerically. Our overall DNS-results are reported in §2.3, while the final asymptotic form of the solution is provided by the ENSresults reported in §2.4. Our asymptotic confirmation of those results is only partially achieved. For though the small-E mainstream analytic solution is

ω = κ 2α k 2 -2α -α2 I 1 (kr/ ) (r/ ) I 1 (k) exp -σ 2 E -1/2 + k 2 τ as τ → ∞ (
(k) I 1 (k) = α (see (2.24b,d
)) characterising (5.1) are only fixed after α is determined by the solution of the E 1/3 -layer problem (see §4), the asymptotic solution of which is possible (see (4.33)) but only useful at extremely small E. For the moderately small E used in the DNS and ENS, it is necessary to solve the shear-layer equations numerically. The §4.2 results determine α = α(E) dependant on E alone and hence complete our hybrid asymptotic-numerical solution.

In contrast to (5.1), the solution (3.8a,b) (cf., (2.23) and (2.25)) for the rigid outer boundary case considered by [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF]

determines ω = 2κ J 1 (j 1 r/ ) (j 1 r/ ) J 2 (j 1 ) exp -σ 2 E -1/2 -j 2 1 τ as τ → ∞ , (5.2)
where j 1 = 3.8317 • • • is a constant, unlike k which is a function of and E being the solution of α(k) = α(E). The derivation of (5.2) is similar to (5.1), particularly (2.23)-(2.27), where now we need Erdélyi et al. 1953, §7.1.4, eq. ( 48)). Furthermore, comparison of the spatial and temporal differences between the solutions (5.1) and (5.2) is striking and reinforces remarks made in §1.1. Though the τ → ∞ limit taken in (5.1) and (5.2) highlights the differences between the stress-free and rigid outer boundary cases, [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] never derived (5.2), because by that late time (even long on the long lateral diffusion time-scale 2 E -1 /Ω) the solution, for all intents and purposes, is negligible. They were concerned with the solution on the shorter spin-down time-scale E -1/2 /Ω on which the solution is still finite. Despite our analysis in §3 of the evolution of the scaled QG-azimuthal angular velocity ω being cast on the 2 E -1 /Ω time-scale, it also pertains to the relatively short time E -1/2 /Ω, provided it is long compared to the time E -1/3 /Ω needed to establish the E 1/3 -layer. It is therefore illuminating to consider the force balances controlling the unscaled angular velocity ω on the shorter spin-down time-scale E -1/2 /Ω by expressing L. Oruba, A. M. Soward and E. Dormy

j1 0 ρ J 1 (ρ) 2 dρ = 1 2 j 2 1 J 2 (j 1 ) 2 (see
(3.1) in its more primitive form

E -1/2 ∂ ω ∂t + σ ω = E 1/2 r 3 ∂ ∂r r 3 ∂ ω ∂r .
(5.3)

From the balance of the Ekman suction term σ ω (bottom friction) on the left-hand side and the lateral diffusion term on the right-hand side, we identify the classical E 1/4 Stewartson (1957) layer. Indeed Barcilon solved (5.3), expressed in boundary layer form (see Barcilon 1968, eqs. (17), ( 18)), subject to Stewartson-like boundary conditions for which there is a steady E 1/4 -layer in the limit E -1/2 t → ∞. [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] also refer to their shear-layer solution on the E -1/2 /Ω time-scale as an E 1/4 -layer. Their primary balance, however, is between the terms on the left-hand side of (5.3) that describe spin-down without identifying any length scale. That fast spin-down is modulated on the slow time τ (see (2.8)), a partition often referred to as two-timing. Then, as the scaled equation (3.1) shows, that modulation diffuses on the lateral diffusion length δ = H(Et) 1/2 , which just happens to equal O(HE 1/4 ) on the spin-down time t = E -1/2 . For t E -1/2 the layer slowly thickens until the solution attains the asymptotic form (5.1). Be that as it may, we reiterate our remarks of §1.2 that the §3.2 power series solution is valid on the spin-down time E -1/2 /Ω and provides our analogue of the rigid boundary solution (Greenspan & Howard 1963, eq. (6.3)).

The cornerstone of our study is the investigation in §4 of the E 1/3 -layer, whose main purpose is the determination of α needed for the QG-mainstream boundary condition dω/dr( ) = αω( ) (see (2.24c)). Nevertheless, amongst the results reported there, the power series solution (4.13a,b) of the zeroth order shear-layer problem in §4.1.2 is of interest in its own right, as it provides an alternative to the usual Fourier series expansion (4.7a,b). Being constructed as the sum of powers of z, each multiplied by functions of the similarity variable Φ = (Ez) -1/3 (r -) (see (4.13), (4.14) and (4.13c)), it is especially accurate in the vicinity of the singularity at (r, z) = ( , 0), where the Fourier series fails or at any rate needs a very large number of terms to converge. Indeed only a very few terms of the power series are needed to give good results throughout most of the shearlayer (see figure 5). It should , however, be stressed that the power series solution is only possible because we have the explicit Fourier forms (4.12a-e) at r = to build on.

Though we relegated our extension of the [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] transient solution of the infinite plane layer problem to appendix A, our new results are both intriguing and useful. The fact that the large time solution (A 5a) involves the factor κ = 1 + O(E 1/2 ) = 1 (see (A 9a)) surprised us, so prompting the question "Where does the O(E 1/2 ) contribution to κ come from?". The possibility, that it originates from the early time behaviour before the Ekman layer has formed, can be dismissed because the Laplace transform solution of appendix A only appeals to a full developed Ekman layer (see (A 2b)). Instead the small O(E 1/2 ) contribution to κ must be a consequence of inertial wave radiation. Certainly inertial waves are visible in the full numerical solution but diminish in size relative to the QG-part as time proceeds. It should be noted that the positive sign of κ -1 ≈ E 1/2 /4 indicates that the QG-flow has decayed less than expected consistent with energy being transferred into inertial waves rather than taken from the QG-flow.

Significantly, the introduction of the factor κ early on in (2.7) has led to remarkable agreement between the transient amplitudes of the DNS-results and the numerical solution of the QG-equation (3.1) displayed in figure 3, which to graph plotting accuracy are virtually indistinguishable except within the E 1/3 -layer. Likewise, the factor µ = 1 -1 2 E 1/2 (see (A 9b)) has proved useful in the definitions (2.32) and (3.9) of f (2q-1) 0 (4τ ) q + α ∞ q=1 g (2q) 0 (4τ ) q (B 7)

F (ζ) = 2 √ π exp( -ζ 2 ), G (
given by (3.15) for n 4 on employing the results below. Substitution of (B 4) with (B 6) into the equations (B 2) and boundary condition (B 3), and then equating coefficients at each respective power of τ leads to At O τ 1/2 , n = 1:

f 10 = 1 2 , g 11 = -1 . ( B 8)
At O τ , n = 2:

f 21 = 1 4 3 2 -α , g 20 = 1 4 3 2 + α , ( B 9a
)

g 22 = - 1 2 3 2 -α . (B 9b)
At O τ 3/2 , n = 3: Though their derivation is straightforward, that for each successive order O τ n/2 becomes increasingly tedious, which is why we have given explicit results for n = 1, • • • , 4.

f 30 =

  .29) similar to the analytic prediction(1.10a,b). The contours of constant ψ DNS and ω DNS for the case = 1, E = 10 -4 , for which q DNS ≈ 96.9648, are illustrated in figure1. A thin Ekman layer is visible near z = 0, linked to the Ekman layer contributions ψ and ω, identified by ψ and ω as τ → ∞. The QG-nature of the mainstream is demonstrated by the ω DNS -contours, which are almost parallel to the z-axis in figure1(a). The mainstream L. Oruba, A. M. Soward and E. Dormy a.b.

Figure 1 .

 1 Figure 1. (Colour online) The contours in the r-z plane of (a) the angular velocity ωDNS and (b) the streamfunction ψDNS, obtained from the DNS of the initial value problem (2.4)-(2.6) for the case = 1, E = 10 -4 in the large τ limit (see (2.29)).

Figure 2 .

 2 Figure 2. Comparison of ENS results ωENS(r) (see (2.32): continuous dark line) and ωENS(r, 1) (continuous light line), normalised by ωENS(0) = ω(0), with the entire analytic mainstream solution ω(r) (dashed line) on 0 r (= 1). Plots are for the cases: (a) E = 3 × 10 -4 ; (b) E = 10 -4 ; (c) E = 3 × 10 -5 ; (d) E = 10 -5 . The values of k, α and σ for each case are given on the respective figures.

Figure 3 .

 3 Figure 3. (Colour online) Profiles of ωDNS(r, τ ) exp -k 2 τ (continuous lines) for the case E = 10 -4 , = 1, for which k 1.9457, α 0.8254, σ 1.0075 (see data on figure 2(b)), at various times τ , together with the corresponding profiles of ω(r, τ ) exp -k 2 τ (dashed lines). The asymptotes ω(0, τ ) exp -k 2 τ ↓ A0 0.7307 and ω( , τ ) exp -k 2 τ ↑ A 1.1355 as τ → ∞ are indicated by the dotted lines. (Inset) The left-hand end-point value log ω(0, τ ) plotted versus τ , together with its large τ asymptote (log A0) + k 2 (log e)τ .

  .12a,b) that solves ∂ω ∂τ = 2 ∂ 2 ω ∂r 2 and meets the initial condition ∂ω ∂r(r, 0) = 0 and the boundary conditions ∂ω ∂r( , τ ) = α and ∂ω ∂r → 0 as ζ → ∞. Integration with respect to r yields ω = 1 + α√ 4τ 1 (ζ) with 1 (ζ) = -erfc ζ (3.13a,b)

Figure 4 .

 4 Figure 4. (Colour online) The right-hand end-point value ω( , τ ) exp -k 2 τ (solid line) plotted versus τ for the case illustrated in figure 3. (a) The dashed curves correspond to ω( , τ ) given by the series expansion (3.15) truncated at various orders: n = 1, Red; 2, Green; 3, Magenta; 4, Blue (see (3.14)). (b) As in (a) but the power series for exp -k 2 τ is employed and the product (3.16b) again approximated at the various orders. We also show the asymptote ω( , τ ) exp -k 2 τ ↑ A 1.1355 as τ → ∞ (dotted line).

  10a,b) equivalently (2.7a,b) and (2.21))

  are simply polynomial expressions. They are rendered unique by the boundary conditions (4.11) and (4.12b,c) at ζ = 0, which require

Figure 5 .Figure 6 .

 56 Figure 5. (Colour online) Contours of (a) streamfunction Ψ and (b) azimuthal velocity V , derived from the complex forms (4.6b,c), in the ζ-z plane. The Fourier series solution (4.7) is red and the similarity solution (4.13) with the power series (4.14) truncated at N = 4 (see (C 1)) is blue. Negative values are identified by dashed lines.

Figure 7 .

 7 Figure 7. (Colour online) As in figure 6 (E = 10 -5 , = 1): Red ENS-contours, (a) ΨENS and (b) VENS (see (4.30a,b)); black ANS-contours, (a) ΨANS and (b) VANS.

Figure 8 .

 8 Figure 8. (Colour online) For the case E = 10 -5 ( ≈ 0.15), = 1, profiles (ANS: black continuous; FS: red dashed; ENS: blue dot-dashed) of (a) 2ΨANS(ζ, 0) = VANS(ζ, 0) and VFS(ζ, 0) versus ζ, and (b) VANS(0, z), VENS(0, z) and VFS(0, z) versus z. (Inset) Blow up (a) near ζ = 0, and (b) near z = 0.

Figure 9 .

 9 Figure 9. (Colour online) As in figure 8, profiles of the scaled z-mean values V = -1 V , (a) VANS(ζ), VENS(ζ) and VFS(ζ), and (b) their derivatives dV /dζ.

  (d). We plot the derivatives dV ANS /dζ(ζ), dV FS /dζ(ζ) = Θ FS (ζ, 0) and dV ENS /dζ(ζ) in figure 9(b).

Figure 10 .

 10 Figure 10. (Colour online) Profile of αANS versus -ln = -(1/6) ln E (black continuous curve).The blue bullets identify the values of αENS ( = 1) for E = 10 -3 , 3 × 10 -4 , E = 10 -4 , 3 × 10 -5 E = 10 -5 respectively (left to right). The best fit power law α = -(1/24π) ln + const. (see (4.33)) is the red dashed line.

  ) evaluated at ζ = const. × for any O(1) constant. However, consideration of the E = 10 -5 case illustrated in figure 9(b) exposes the folly of this expectation at finite E. In the vicinity of ζ = --0.15, the value of dV FS /dζ varies rapidly and the choice of different O(1) constants will determine enormously different values of α FS . Nevertheless as a somewhat futile gesture, we identified in figure 10 the power law α = -(4π) -1 ln E + const. (4.33)

  ζ) = erfc ζ = express each n (ζ) in the solution (3.14) for ω in the formn = f n (ζ)F (ζ) + g n (ζ)G(ζ) , n/2 f n (ζ)F (ζ) + g n (ζ)G(ζ) , (B 5)and seek polynomial solutionsf 1 (ζ) = f 10 , g 1 (ζ) = g 11 ζ , (B 6a) f 2 (ζ) = f 21 ζ , g 2 (ζ) = g 20 + g 22 ζ 2 , (B 6b)f 3 (ζ) = f 30 + y 32 ζ 2 , g 3 (ζ) = g 31 ζ + g 33 ζ 3 , (B 6c) f 4 (ζ) = f 41 ζ + f 43 ζ 3 , g 4 (ζ) = g 40 + g 42 ζ 2 + g 44 ζ 4 , (B 6d )in which f nm and g nm are constants. The polynomials for n > 4 can be expressed in the obvious way. They determine the end point (r = ) value ω( , τ ) = 1 + α √ πτ ∞ q=1

† Email address for correspondence: andrew.soward@ncl.ac.uk

Acknowledgements

This project was partially supported by the French National Center for Scientific Research (CNRS) under the interdisciplinary grant Inphyniti/MI/CNRS. The numerical simulations were performed using HPC resources from GENCI-IDRIS (Grant 2015-100584 and 2016-100610).

Appendix A. Final decay of QG-motion in an unbounded plane layer

We extract the salient details of solution of the spin-up problem considered by [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] for two parallel coaxial infinite disks in their section 3 and extend them as needed for the infinite plane layer solution of our initial value spin-down problem (2.4-2.6). Like [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] we assume the similarity form

and solve by the Laplace transform

method. We define Greenspan & Howard 1963, eq. (3.6), in which (m 1 , m 2 ) → (m + , m -)) so that, provided that |m ± | = O E -1/2 1, the solution may be written

where Greenspan & Howard 1963, particularly eq. (3.11)).

The E -1/2 z → ∞ limit of (A 2b) determines the mainstream values

The QG-flow for t = O E -1/2 1 is determined by the residue at the pole p = p 0 , where D(p) = 0, of the inverse Laplace transform. It yields

where Assuming that p = O E 1/2 , we obtain

We set p 0 = -E 1/2 σ (real σ) and substitute into p

with the solution expansion

In addition (A 5b), (A 6b) and (A 7a,b) determine

Appendix B. The small τ power series solution

To solve (3.1)-(3.4) for τ 1, we make the preliminary change of variables

where Appendix C. Power series for Ψ and V for N 4

We provide explicit series solutions for Ψ = Im{Ξ} and V = Re{Υ }, where Ξ and Υ are defined by (4.13a,b). We truncate the series (4.14a,b) for Ξ N (Φ, z), Υ N (Φ, z) that approximate Ξ(Φ, z), Υ (Φ, z) at N = 4 to obtain