The propagation and decay of a coastal vortex on a shelf
TThis draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics The propagation and decay of a coastalvortex on a shelf
Matthew N. Crowe and Edward R. Johnson
Department of Mathematics, University College London, London, WC1E 6BT, UK(Received xx; revised xx; accepted xx)
A coastal eddy is modelled as a barotropic vortex propagating along a coastal shelf. Ifthe vortex speed matches the phase speed of any coastal trapped shelf wave modes, ashelf wave wake is generated leading to a flux of energy from the vortex into the wavefield. Using a simply shelf geometry, we determine analytic expressions for the wave wakeand the leading order flux of wave energy. By considering the balance of energy betweenthe vortex and wave field, this energy flux is then used to make analytic predictionsfor the evolution of the vortex speed and radius under the assumption that the vortexstructure remains self similar. These predictions are examined in the asymptotic limit ofsmall rotation rate and shelf slope and tested against numerical simulations.If the vortex speed does not match the phase speed of any shelf wave, steady vortexsolutions are expected to exist. We present a numerical approach for finding thesenonlinear solutions and examine the parameter dependence of their structure.
Key words:
Waves in rotating fluids, Topographic effects, Vortex dynamics
1. Introduction
The interaction of interior ocean flows with coastal boundaries is a complicated multi-scale problem with important implications for the dissipation of mesoscale energy and thegeneration of potential vorticity (Dewar et al. et al. et al. et al. .2006; Dewar & Hogg 2010;Hogg et al. et al. et al. a r X i v : . [ phy s i c s . f l u - dyn ] F e b M. N. Crowe & E. R. Johnson travelling vortices would similarly generate waves, however, since the only source ofenergy for the wave-field is the kinetic energy of the vortex, the formation of a wave-field would lead to a loss of vortex energy and hence a decay of the vortex. This leadsto a feedback mechanism where the generated modes depend on the properties of thevortex and the vortex decay depends on the wave energy flux. Flierl & Haines (1994)used an adjoint method to examine the decay of a modon on a beta plane. They foundthat as the vortex decayed, mass was ejected from the rear. Therefore, unlike energy,momentum and enstrophy were not conserved between the vortex and wave-field. Thevalue of the maximum vorticity was argued to be a second conserved quantity and usedto make analytical predictions for the decay of the modon speed and radius. Johnson& Crowe (2021); Crowe et al. (2021) estimated the decay of the Lamb-Chaplygin dipole(Meleshko & van Heijst 1994) and Hill’s vortex (Hill 1894) in rotating and stratified flowsby calculating the work done by the leading order wave drag and equating this to the lossof vortex energy. Conservation of maximum vorticity was again used to close the systemand shown to be valid using numerical simulations.Here we consider a simple analytical model of a moving vortex on the boundary of acoastal shelf with the aim of determining the long term evolution. Our vortex is takento consist of a near semi-circular region of vorticity centred on the boundary. Using themethod of images, this may be modelled as a dipolar vortex with the dipole strengthdetermined by the vortex speed and radius. We begin in Section 3 by considering thegeneration of shelf waves by a moving vortices. As expected, we observe that a wave-fieldwill only be generated if the vortex speed matches the phase speed of any shelf wavesand hence vortices moving faster than, or in the opposite direction to, every shelf wavemode will not generate a wave wake. The generation of these waves leads to a flux ofenergy from the vortex to the wave-field resulting in a decay of the vortex. We use asimple energy balance to estimate this decay and present analytical results for the caseof asymptotically small rotation rate and shelf slope. Our predictions are tested againstnumerical simulations in Section 4. Finally, in Section 5, we examine the case where thevortex does not generate waves. We expect steady vortex solutions to exist and presenta numerical approach for finding these fully nonlinear solutions.
2. Setup
Our starting point is the two-dimensional rotating shallow water equations under therigid lid assumption. Let Ox (cid:48) y be Cartesian coordinates, fixed in the topography, where x (cid:48) describes the distance along a straight coastline and y describes the distance in theoffshore direction. We consider a near-semicircular vortex moving along the coastline withspeed U ( t ) and introduce coordinates following the vortex centre (Johnson & Crowe 2021)by defining x = x (cid:48) − (cid:90) t U ( t (cid:48) ) d t (cid:48) , (2.1)and working in the coordinate system Oxy . All quantities are specified by their valuein the frame of the topography unless stated otherwise. We thus simply work withtopographic frame variables expressed as functions of the moving coordinates ( x, y ) ratherthan working in variables relative to a frame translating with speed U ( t ). The advantageof this formulation is that fictitious forces, that would arise from treating quantitiesrelative to the accelerating vortex frame, are absent. Throughout, whenever consideringvortices which are semi-circular, or asymptotically close to semi-circular, we will denotethe radius by a ( t ). ropagation and decay of a shelf vortex U(t) a(t) x′ y H(y)D ε = / Ro x (a) Figure 1: Our non-dimensional setup showing a vortex of radius a ( t ) moving along acoastal boundary with speed U ( t ). The layer depth, H ( y ), is shown as a shelf region ofincreasing depth of width D matched to a constant depth ocean. The system is rotatingwith inverse Rossby number of (cid:15) .The non-dimensional equations in terms of ( x, y ) governing the horizontal velocitycomponents ( u, v ) relative to the topography are thus ∂u∂t − U ∂u∂x + u ∂u∂x + v ∂u∂y − (cid:15)v = − ∂p∂x , (2.2 a ) ∂v∂t − U ∂v∂x + u ∂v∂x + v ∂v∂y + (cid:15)u = − ∂p∂y , (2.2 b ) ∂∂x ( uH ) + ∂∂y ( vH ) = 0 , (2.2 c )where (cid:15) = 1 /Ro is the inverse Rossby number and H is the layer depth. We take H torepresent a shelf with the depth varying only in the offshore, y , direction so H = H ( y ).Our setup is shown in Fig. 1 for a vortex of radius a ( t ) and a depth profile, H ( y ),consisting of a shelf of width D joined to a region of constant depth.Eq. (2.2) can be combined to give a single evolution equation for the potential vorticity (cid:18) ∂∂t − U ∂∂x (cid:19) q + 1 H J [ ψ, q ] = 0 , (2.3)where the velocity can be expressed using a volume flux streamfunction( u, v ) = 1 H (cid:18) − ∂ψ∂y , ∂ψ∂x (cid:19) , (2.4)the vorticity is given by ζ = ∂v∂x − ∂u∂y = 1 H ∂ ψ∂x + ∂∂y (cid:20) H ∂ψ∂y (cid:21) , (2.5)and q denotes the potential vorticity (PV) in the layer q = ζ + (cid:15)H . (2.6) M. N. Crowe & E. R. Johnson
We impose a wall at y = 0 with the boundary condition v ( x,
0) = 0 or equivalently ψ ( x,
0) = 0. Far from the wall, disturbances are assumed to decay so ( u, v ) → y → ∞ .Throughout this study we will consider a shelf profile for H consisting of a shelf regionwith exponentially decreasing depth matched to a constant depth ocean. Where possible,general results will also be given in terms of the arbitrary profile H = H ( y ). Our shelfprofile is given by H ( y ) = (cid:40) exp[ βy ] , y (cid:54) D, exp[ βD ] , y (cid:62) D, (2.7)where D describes the shelf width and β describes the shelf slope.We take our vortex solution to be dipolar and centred at ( x, y ) = (0 ,
0) in our vortex-following coordinates. This corresponds to a single vortex in y > µ directed in the positive x direction.Therefore, far from the vortex we have ψ = µ δ ( x ) at y = 0 + . (2.8)In the limit of small (cid:15) and β , our vortex solution is given by the classical Lamb-Chaplygindipole (Meleshko & van Heijst 1994) and hence ψ = (cid:40) − U y + ( Kr ) y J ( Ka ) Kr , r < a − U a y/r , r > a, (2.9)for r = x + y . Here a is the semi-circular vortex radius, J and J are Bessel functionsof the first kind and K = j /a where j ≈ . . Inthis case the dipole strength may be calculated as µ = 2 πU a . A numerical approach forfinding these vortex solutions for a general depth profile, H , and rotation rate, (cid:15) , is givenin Section 5.
3. Vortex decay
This shelf system admits shelf wave solutions so if the vortex is travelling in the samedirection as the phase speed of these waves (requiring (cid:15)U > et al.
The topographic wave field
We begin by determining the linear topographic wave solutions admitted by the systemin the absence of a vortex. Our wave equation is obtained by linearising Eq. (2.3) andsetting U = 0 to get ∂∂t (cid:20) ∇ ψ − H y H ∂ψ∂y (cid:21) − (cid:15)H y H ∂ψ∂x = 0 . (3.1) ropagation and decay of a shelf vortex H to be our shelf profile from Eq. (2.7) and assume wavelike solutions ofthe form ψ = √ H ˜ φ ( y ) exp(i ωt − i kx ) to obtain (cid:20) ∂ ∂y − k (cid:21) ˜ φ = (cid:40) − κ ˜ φ y (cid:54) D, y (cid:62) D, (3.2)where κ ( k, ω ) = (cid:15)βkω − β . (3.3)We take boundary conditions of ˜ φ = 0 on y = 0, ˜ φ → y → ∞ and ˜ φ , ˜ φ y arecontinuous across y = D . For y > D our solution is of the form˜ φ = C ( k ) exp ( −| k | y ) , (3.4)for some C hence we can impose the boundary condition˜ φ y + | k | ˜ φ = 0 on y = D, (3.5)and only consider the shelf region y ∈ [0 , D ]. In the shelf region we have solution givenby ˜ φ = C ( k ) sin (cid:104)(cid:112) κ − k y (cid:105) , (3.6)where the square root term make be complex. Using the boundary condition at y = D ,Eq. (3.5), we obtain the dispersion relationtan (cid:104)(cid:112) κ − k D (cid:105) = − √ κ − k | k | , (3.7)with solution describing a countably infinite set of modes with differing wave numberand offshore structure. To proceed we define l = (cid:112) κ − k , (3.8)so l can be thought of as the offshore wavenumber discretised by the shelf boundary at y = D . We can now solve Eq. (3.7) numerically for l ( k ) for each mode and plot thefrequency, ω , and phase speed, c p = ω/k , by combining Eq. (3.3) and Eq. (3.8) to get ω = (cid:15)βkk + l + β / c p = (cid:15)βk + l + β / . (3.9)For (cid:15) >
0, the phase speed of the waves is positive for all wavenumbers and conversely for (cid:15) <
0, the topographic wave phase speed is always negative. Additionally, the frequencyis odd in k , so we only need to consider waves with k (cid:62)
0. We note that the offshorewavenumber for a given mode, l ( k ), is an increasing function of k and lies within theinterval l ∈ (cid:18)(cid:20) n − (cid:21) πD , nπD (cid:19) , (3.10)for mode number n = { , , , . . . } with l (0) = (cid:20) n − (cid:21) πD and l ( k ) → nπD as k → ∞ . (3.11)Since ω = c p k , the group velocity is given by c g = ∂ω∂k = c p + k ∂c p ∂k , (3.12) M. N. Crowe & E. R. Johnson k n = 1n = 2n = 3n = 4n= 5 (a) k c p n = 1n = 2n = 3n = 4n= 5 (b) Figure 2: Plots of the frequency, ω , (a) and phase speed, c p , (b) for the first five modeswith (cid:15) = 0 . β = 0 . D = 25 . k (cid:54) = 0 hence c g < c p for k (cid:54) = 0.Fig. 2 shows the frequency and phase speed for the first five modes with (cid:15) = 0 . β = 0 . D = 25 .
6. These curves are consistent with classical results for topographicRossby waves. For a vortex to generate a wave field, the speed of the vortex must matchthe phase speed of one or more waves hence a vortex cannot generate any waves if itmoves faster than the fastest mode or moves in the opposite direction to the topographicwaves. Therefore, for our choice of topography, a vortex moving with speed U will onlygenerate waves if 0 < (cid:15)U < (cid:15) βl (0) + β / , (3.13)where l (0) = π/ (2 D ) is the offshore wavenumber of the lowest mode for k = 0. We havemultiplied this condition through by (cid:15) to ensure it holds for both positive and negative (cid:15) .We note that c g (cid:54) c p for all modes hence for a vortex moving with speed U = c p , energywill be emitted from the rear of the vortex. This radiation condition will be requiredlater.If the vortex does not generate waves we expect that steady vortex solutions will exist,these solutions can be found using the method outlined in Section 5.3.2. Waves generated by a moving vortex
We now determine the amplitude of the vortex generated wavefield. Working incoordinates following the vortex ( U (cid:54) = 0) and looking for steady wave solutions we obtainthe linearised wave equation − U ∂∂x (cid:20) ζ + (cid:15)H (cid:21) − (cid:15)H y H ∂ψ∂x = 0 . (3.14)Substituting for ψ = √ Hφ gives ∇ φ = (cid:40)(cid:16) β − (cid:15)βU (cid:17) φ y (cid:54) D, y (cid:62) D, (3.15)where we note that κ ( k, U k ) = (cid:15)βU − β . (3.16) ropagation and decay of a shelf vortex κ must be positive if any waves are generated by the vortex by Eq. (3.13). Farfrom the vortex, the vortex appears as a point dipole of strength µ hence we impose theboundary condition φ = µ δ ( x ) at y = 0 + . (3.17)We also take φ → y → ∞ and impose continuity of φ and φ y across y = D .We now express φ using a Fourier transform as φ ( x, y ) = 12 π (cid:90) ∞−∞ (cid:98) φ ( k, y ) exp(i kx ) d k, (3.18)where (cid:98) φ satisfies the system (cid:20) ∂ ∂y − k (cid:21) (cid:98) φ = (cid:40) − κ (cid:98) φ y (cid:54) D, y (cid:62) D, (3.19)subject to (cid:98) φ = µ y = 0 , (3.20) (cid:98) φ → y → ∞ , (3.21)and the requirement that (cid:98) φ and (cid:98) φ y are continuous across y = D . The solution for (cid:98) φ isgiven by (cid:98) φ = µ (cid:16) (cid:98) C sin (cid:2) √ κ − k y (cid:3) + cos (cid:2) √ κ − k y (cid:3)(cid:17) y (cid:54) D, (cid:16) (cid:98) C sin (cid:2) √ κ − k D (cid:3) + cos (cid:2) √ κ − k D (cid:3)(cid:17) e −| k | ( y − D ) y (cid:62) D, (3.22)where (cid:98) C = √ κ − k sin (cid:2) √ κ − k D (cid:3) − | k | cos (cid:2) √ κ − k D (cid:3) √ κ − k cos (cid:2) √ κ − k D (cid:3) + | k | sin (cid:2) √ κ − k D (cid:3) . (3.23)We note that the denominator of (cid:98) C vanishes if the dispersion relation in Eq. (3.7) issatisfied therefore the wave modes correspond to the residues of these poles.By the radiation condition that c g − U < x < | x | (cid:29)
1. For large x , the exponential term in Eq. (3.18) is strongly oscillatoryand the terms without poles decay as 1 /x . We therefore consider only the terms containing (cid:98) C and form a closed contour by including the arc | k | = R with R → ∞ in the lower halfof the complex k plane. All poles occur along the real line and are taken to lie within thecontour. Finally, the contribution from the arc vanishes as R → ∞ giving that φ ∼ − i µ (cid:88) Res[ (cid:98) φ, k n ] exp(i k n x ) , (3.24)where Res[ f, x ] denotes the residue of f and x . Here we sum over all poles of (cid:98) φ and havegained an additional factor of − (cid:98) φ is even in k there will be a pole at − k = k n for each pole at k = k n withresidue of the opposite sign. We therefore have φ ∼ i µ N (cid:88) n =1 Res[ (cid:98) φ, k n ] [exp( − i k n x ) − exp(i k n x )] = µ N (cid:88) n =1 Res[ (cid:98) φ, k n ] sin( k n x ) , (3.25)where the k n are the positive poles of (cid:98) C and hence are the solutions to the dispersion M. N. Crowe & E. R. Johnson k c p n = 1n = 2n = 3n = 4n= 5 (a) k c p n = 1n = 2n = 3n = 4n= 5 (b) Figure 3: Plots of the phase speed, c p , for the first five modes with β = 0 . D = 25 . (cid:15) = 0 .
2. (b) (cid:15) = 1. The dotted line denotes U = 1 and the k n are determined as theintersections c p = U and shown by the black circles. We observe one mode for (cid:15) = 0 . (cid:15) = 1 for this choice of parameters.relation corresponding to a mode of phase speed U . N describes the number of modesand will be determined later. By differentiating the denominator of (cid:98) C we find that theresidues are given byRes[ (cid:98) φ, k n ] = √ κ − k n k n D sin (cid:104)(cid:112) κ − k n y (cid:105) y (cid:54) D, √ κ − k n k n D sin (cid:104)(cid:112) κ − k n D (cid:105) e − k n ( y − D ) y (cid:62) D, (3.26)and hence φ ∼ µ (cid:80) Nn =1 (cid:20) √ κ − k n k n D sin (cid:104)(cid:112) κ − k n y (cid:105) sin( k n x ) (cid:21) y (cid:54) D, (cid:80) Nn =1 (cid:20) √ κ − k n k n D sin (cid:104)(cid:112) κ − k n D (cid:105) sin( k n x ) e − k n ( y − D ) (cid:21) y (cid:62) D. (3.27)This result may alternatively be written in terms of an offshore wavenumber l n = (cid:112) κ − k n satisfying tan( l n D ) = − l n k n , (3.28)as φ ∼ µ (cid:80) Nn =1 (cid:104) l n k n D sin( l n y ) sin( k n x ) (cid:105) y (cid:54) D, (cid:80) Nn =1 (cid:104) l n k n D sin( l n D ) sin( k n x ) e − k n ( y − D ) (cid:105) y (cid:62) D. (3.29)The values of ( k n , l n ) can be determined numerically by finding the roots of c p = U usingEq. (3.9) or solving Eq. (3.28) directly. Finally, by calculating the maximum value of c p for each mode and comparing this to U , we may determine the total number modes as N = (cid:36)
12 + Dπ (cid:114) (cid:15)βU − β (cid:37) . (3.30)The case of N = 0 corresponds to the vortex moving faster than the fastest wave and isequivalent to the second inequality in Eq. (3.13) not being satisfied. ropagation and decay of a shelf vortex k n , of the modes, determined by solving c p = U for β = 0 . D = 25 . U = 1 and (cid:15) = 0 . ,
1. The value of l n can now be easilydetermined using l n = (cid:112) κ − k n . We note that if a vortex slows down, it will generatean increased number of modes.3.3. Wave energy flux and vortex decay
As the vortex generates waves, it loses energy to the wavefield and therefore decays.This energy flux is given to leading order by the quadratic pressure work and the totalenergy flux can be determined by calculating the flux across the line x = − L for L (cid:29) F = (cid:90) ∞ Hpu d y, (3.31)where the factor of H is obtained by integrating over the layer depth. For our linear wavefield the pressure is given by p = U u to leading order and hence F = U (cid:90) ∞ Hu d y = U (cid:90) ∞ H (cid:18) ∂ψ∂y (cid:19) d y = U (cid:90) ∞ φ y + H y H φφ y + H y H φ d y. (3.32)We may now calculate F using Eq. (3.29). We note that since φ depends on x , ourtotal energy flux will depend also on x , however we can average over x to obtain thetime-averaged energy flux. This is equivalent to time-averaging the energy flux over thetimescale of the wave generation which is assumed to be much shorter than the timescaleof the vortex decay. Denoting this average by (cid:104)∗(cid:105) we use that (cid:104) sin( ax ) sin( bx ) (cid:105) = (cid:40) a (cid:54) = b a = b, (3.33)so quadratic quantities in φ may be written as a sum of the terms of φ squared similarlyto an orthogonality approach.Our final energy flux results is F = U µ N (cid:88) n =1 l n (1 + k n D ) (cid:20) D (cid:18) (cid:15)βU − k n (cid:19) + β (cid:18) βk n l n + 1 (cid:19) sin ( l n D ) (cid:21) , (3.34)which we note is always positive. Equating this flux with the loss of vortex energy, E ,gives dEdt = − F. (3.35)If we now assume that the vortex remains self similar throughout the evolution, theenergy, E , and dipole strength, µ , may be determined in terms of the vortex speed U and radius a . We now have two quantities, U and a , with a single evolution equation,Eq. (3.35), so a second equation is required to close the system. Following Flierl & Haines(1994), Johnson & Crowe (2021) and Crowe et al. (2021) we choose conservation of centrevorticity so that the maximum vorticity of the vortex remains constant throughout theevolution. This gives a second equation dη c dt = 0 , (3.36)where η c = η c ( U, a ) is the maximum vorticity within the vortex can be determined fromthe vortex solution, either numerically or analytically.0
M. N. Crowe & E. R. Johnson U FF b (a) U FF b (b) Figure 4: Plots of the wave energy flux, F , as a function of U for β = 0 . D = 25 . µ = 2 πU and (cid:15) = 0 . (cid:15) = 1 (b). If the vortex speed, U , exceeds the fastest wavethere is no wave energy flux, this occurs for U > .
19 for (cid:15) = 0 . U > .
96 for (cid:15) = 1. The red line shows the predicted flux bound, F b , where F ≈ F b if the alongshorewavenumber, k n , vanishes.In the case of asymptotically small β and (cid:15) the vortex solution reduces the the classicalLamb-Chaplygin dipole and the quantities may be determined analytically from Eq. (2.9)as E ( U, a ) = πU a , µ ( U, a ) = 2 πU a , and η c ( U, a ) ∝ Ua . (3.37)Therefore Eqs. (3.35) and (3.36) give ddt ( U a ) = − πU a N (cid:88) n =1 l n (1 + k n D ) (cid:20) D (cid:18) (cid:15)βU − k n (cid:19) + β (cid:18) βk n l n + 1 (cid:19) sin ( l n D ) (cid:21) , (3.38)and ddt (cid:18) Ua (cid:19) = 0 . (3.39)Therefore a ( t ) ∝ U ( t ) so Eq. (3.38) may be solved as an equation for U ( t ) subject tosome initial condition ( U, a ) = ( U , a ) at t = t . (3.40)We note that since the wavevector, ( k n , l n ), and number of modes, N , both have acomplicated dependence on U this equation would have to be solved numerically.Fig. 4 shows the wave energy flux, F , as a function of U for β = 0 . D = 25 . a = 1, µ = 2 πU a and (cid:15) ∈ { . , } . For large values of U there is no wave energy flux as thereare no waves which match the vortex speed. As we decrease U , more waves can matchthe vortex speed resulting in a non-zero wave energy flux. Each new wave mode appearswith k n = 0 resulting in the spikes shown in Fig. 4. When k n = 0 we have l n = κ andthe wave energy flux of this new mode dominates the energy flux expression giving F ∼ F b = U µ (cid:18) (cid:15)βU − β (cid:19) (cid:32) (cid:15)βDU + β sin (cid:34)(cid:114) (cid:15)βU − β D (cid:35)(cid:33) . (3.41)As the vortex speed is reduced further, the alongshore wavenumber of the newly generatedmode increases and the associated energy flux decreases as the system moves away from ropagation and decay of a shelf vortex The small U limit Assuming that for a given speed U , the energy flux is dominated by the contributionfrom the longest alongshore mode (with smallest k n and largest L n ), we have F ≈ U µ l N (1 + k N D ) (cid:20) D (cid:18) (cid:15)βU − k N (cid:19) + β (cid:18) βk N l N + 1 (cid:19) sin ( l N D ) (cid:21) , (3.42)since k N < k N − < · · · < k . We now take U to be small and consider 2 casescorresponding to different ratios between the alongshore wavelength, λ n = 2 π/k n , andthe shelf width, D .Firstly we consider the case when the longest alongshore mode has just appeared. Herewe have k N D (cid:28) l N ∼ (cid:112) (cid:15)β/U so the energy flux may be estimated as F ∼ µ (cid:15) β DU , (3.43)so taking (cid:15) , β small we may use Eq. (3.37) to determine approximate solutions for U ( t )and a ( t ) as U ( t ) U , a ( t ) a ∼ (cid:18) (cid:15) β Dπ a U ( t − t ) (cid:19) − . (3.44)Conversely, if the longest alongshore mode has both k N and l N large, ( k N , l N ) ∼ U − / ,we have k N D (cid:29) F ∼ C (cid:15)βµ D . (3.45)Here C is a weakly time dependent quantity determined by k N U / and l N U / . Taking (cid:15) and β to be small we may again use Eq. (3.37) to obtain U ( t ) U , a ( t ) a ∼ (cid:18) (cid:15)βC π D a U ( t − t ) (cid:19) − / , (3.46)where we’ve assumed C is approximately constant close to some initial time, t .These solutions therefore correspond to (piecewise) polynomial decay of speed andradius similar to the case of the Rossby wave decay of a beta plane modon consideredby Flierl & Haines (1994) and Johnson & Crowe (2021). We note that the exponent islarger for the case where the longest alongshore mode has a very small value of k N whichis consistent with our observations from Fig. 4 that the energy flux is largest just after anew mode has appeared.
4. Numerical simulations
To test our predictions we perform numerical simulations using Dedalus (Burns et al. U f , in the along-shore( x ) direction. We use the numerical domain ( x, y ) ∈ [ − . , . × [0 , .
2] with 1024gridpoints in each direction and decompose fields in terms of a Fourier basis in the x direction and a compound Chebyshev basis in the y direction with separate Chebyshevexpansions on and off the shelf. Solutions are integrated for t ∈ [0 ,
50] using a secondorder semi-implicit BDF scheme with a timestep of 10 − . We take boundary conditions2 M. N. Crowe & E. R. Johnson of no flow through the walls at y = 0 and y = 51 . ν = 1 . × − for numerical stability. We note that the use of a Fourierbasis in the x direction results in a periodic boundary and hence waves may loop aroundthe domain and interfere with the vortex. However, the stop time, t = 50, is found tobe sufficiently early that these waves do not interact with the vortex. Similarly, the solidwall at y = 51 . y > D ),we’d expect any effects of this rigid wall to be exponentially small.For all simulations the shelf slope, β , and shelf width, D , are chosen as ( β, D ) =(0 . , . | U (0) | = 1 and radius a (0) = 1.The frame speed, U f , is set to match the speed of this initial vortex. Therefore, we expectthe vortex to remain close to x = 0 throughout the evolution with deviations occurringas the vortex speed changes.The effects of non-zero β and (cid:15) are to modify the initial vortex leading to a transientadjustment phase at the beginning of the simulation as the vortex adjusts to the effectsof rotation and shelf slope and the wave field begins to develop. For small (cid:15) and β , thisadjustment is small and the vortex remains approximately a Lamb-Chaplygin dipole witha modified speed and radius. In order to compare with our theoretical predictions, wetake the values of U and a to be the speed and radius after this adjustment phasewith t = t describing the time taken for this adjustment to occur. A value of t = 1 isfound to be sufficient and comparison is made with the theory for t (cid:62) t . We note thataccurately determining the values of U and a from the numerical data is difficult. Thiscan present issues when comparing with our theoretical predictions due to the sensitivedependence of Eq. (3.38) on these quantities. The vortex speed, U , is determined bytracking the position of the vorticity maximum and the vortex radius, a , is estimatedusing the point at which the vorticity becomes 2% of its maximum value.Fig. 5 shows the streamfunction, ψ , for the final timestep, t = 50, of our numericalsimulations for a range of parameters. Panels (a) and (b) show vortices which arerespectively moving faster than and in the opposite direction to all shelf wave modes.For these simulations, the value of ψ c is conserved to within the error expected due toviscous effects and while a very weak wave signature is observed, this is likely the resultof transient waves generated during the initial adjustment. Fig. 5.(c) shows ψ ( x, y, U (0) = 1 and (cid:15) = 0 . λ n = 2 π/k n = 66 .
4. This wavenumber approximately matchesthe observed wave which we note is likely to be restricted by the length of the domain.Finally Fig. 5.(d) shows ψ ( x, y,
50) for U (0) = 1 and (cid:15) = 1. The initial speed, U (0) = 1,matches the phase speed of the first three modes however we observe that the vortexundergoes significant adjustment due to the fairly large value of (cid:15) and adjusts to avalue of U ≈ .
1. This value of U only matches the speed of the first two modesand, according to our theoretical predictions, corresponds to modes with wavelengthsof λ n = 22 . λ n = 30 . l n = 0 .
11 and l n = 0 . (cid:15) and β we can show using the Lamb-Chaplygin solution that the totalstreamfunction at the position of the maximum vorticity is proportional to U a . Therefore, ropagation and decay of a shelf vortex (a) (b)(c) (d) Figure 5: The streamfunction, ψ , as a function of position in a frame moving with thespeed of the initial vortex, U f = U (0). Results are shown for β = 0 . D = 25 . t = 50, for various inverse Rossby numbers, (cid:15) , and initial vortex speeds, U (0). The vortexadjusts slightly due to finite (cid:15) and β effects so the value of U = U taken at t = t = 1can differ from U (0) by up to 10 − (cid:15), U (0)) = (0 . , (cid:15), U (0)) = (0 . , − (cid:15), U (0)) = (0 . , (cid:15), U (0)) = (1 , ψ c ( t ) ψ = U ( t ) a ( t ) U a = U ( t ) U , (4.1)where ψ c is the value of the streamfunction at the position of maximum vorticity and ψ is the value of ψ c at t = t = 1. The normalised value of ψ c can be easily determinedfrom our simulations and used to test our decay predictions by comparing with solutionsof Eqs. (3.38) and (3.39).For very small values of (cid:15) and β the wave energy flux, F , is small such that the vortexdecay, and hence the decrease in ψ c , is slow. Since the effect of viscosity is to decreasethe domain averaged energy by around 1 −
2% over the time interval t ∈ [0 , M. N. Crowe & E. R. Johnson (a) (b)
Figure 6: Plots of the streamfunction, ψ , showing the formation of the wavefield for (cid:15) = 0 . U = 1 . a = 1, β = 0 . D = 25 . t , t = 25 (a) and t = 50(b). Similarly to Fig. 5, the solution is shown in a frame moving with speed U f = 1 inthe x direction. t c TheorySimulation (a) t c (b) Figure 7: (a) The normalised value of ψ c as a function of t , the solid line showsour numerical results and the dashed line gives our analytical prediction from solvingEqs. (3.38) and (3.39). (b) The normalised value of η c from our numerical simulation asa function of time. All results are shown for (cid:15) = 0 . U = 1 . a = 1, β = 0 . D = 25 . (cid:15) greater that 1, while the wavefield remains small due to small β ,the vortex is no longer well described by the Lamb-Chaplygin solution and the asymptoticexpressions for E , µ and η c in Eq. (3.37) begin to deviate from the true values. Thoughthese deviations are fairly small despite an order 1 value of (cid:15) , we find that Eq. (3.35)is very sensitive to the values of E and µ and our prediction gives only the order ofmagnitude of the decay scale rather than an accurate result. As a compromise betweenthese limits we consider here the case of (cid:15) = 0 . β = 0 . U (0) , a (0)) = (1 ,
1) whichwe observe is well described by the Lamb-Chaplygin solution. After the initial adjustmentwe find that U = 1 . a = 1 corresponding to two shelf wave modes with alongshorewavelengths of λ n = 25 . λ n = 38 . ropagation and decay of a shelf vortex ψ , from our numerical simulation with ( (cid:15), β, D, U , a ) =(0 . , . , . , . , k n and smaller offshore wavenumber l n ), propagates slower in the x direction due to amore negative value of the phase speed, c g − U f , and can be seen looping around thedomain ahead of the slower mode. Fig. 7 shows the values of ψ c and η c from the samesimulation. In Fig. 7.(a) we observe fairly close agreement between our numerical valueof ψ c and our theoretical prediction using Eqs. (4.1) and (3.38). The numerical value of ψ c appears to oscillate around our theoretical prediction; this is likely a consequence ofthe time-averaging of the wave energy flux in Eq. (3.32) which smooths the theoreticalprediction. In Fig. 7.(b) we plot the value of the maximum vortex vorticity, η c , as afunction of time in order to test our assumption (see Eq. (3.36)) that this quantity isconserved over the decay scale of the vortex. We observe that η c decreases by around2% over the course of the simulation and since this is much smaller than the decreasein ψ c and similar in magnitude to the effects of viscous dissipation we believe that thisassumption is valid.A supplementary movie file (Movie.mp4) is included to show the temporal evolutionof the wave field from the simulation with ( (cid:15), β, D, U , a ) = (0 . , . , . , . ,
1) shownin Figs. 6 and 7. We plot the mass fluxes, (
Hu, Hv ), as functions of time, t , and position,( x, y ), in a frame moving with constant speed, U f = 1. The formation of two modeswith different structure and speed can be clearly seen. Additionally, we observe that thewaves looping around the domain due to the periodic x direction are unlikely to have asignificant effect on the vortex for t (cid:46)
50. Further, the use of rigid boundary at y = 51 . y = 51 .
5. The steady nonlinear problem
We have shown that if the speed of the vortex matches the topographic wave speedfor some wavenumber then the vortex will generate waves and decay due to the transferof energy from the vortex to the wave field. If however there is no wave speed whichmatches the vortex speed, we expect steady vortex solutions to exist. This can occur ifthe vortex moves in the opposite direction to the waves or moves faster than the fastesttopographic wave. The conditions for a decaying vortex are given in Eq. (3.13) and wewill focus here on how to determine steady vortex solutions to the full nonlinear systemwhen these conditions are not satisfied.Neglecting the time derivative in Eq. (2.3) and combining the constant advection termwith the Jacobian gives J (cid:20) ψ + U (cid:90) H d y, ζ + (cid:15)H (cid:21) = 0 , (5.1)hence we have that the potential vorticity can be written as a function of the totalstreamfunction as ζ + (cid:15)H = F (cid:18) ψ + U (cid:90) H d y (cid:19) . (5.2)The function F may now be determined outside the vortex using the far field conditionthat ζ, ψ → x → ∞ and hence F (cid:18) U (cid:90) H ( y ) d y (cid:19) = (cid:15)H ( y ) , (5.3)6 M. N. Crowe & E. R. Johnson for all y . Defining A ( y ) = (cid:90) y H ( y (cid:48) ) d y (cid:48) , (5.4)as the cross-sectional area in the offshore region [0 , y ] gives that F ( z ) = (cid:15)H ( A − ( z/U )) . (5.5)The full nonlinear problem outside the vortex can now be written as1 H ∂ ψ∂x + ∂∂y (cid:20) H ∂ψ∂y (cid:21) + (cid:15) = (cid:15)H ( y ) H ( A − ( ψ/U + A ( y ))) , (5.6)and solved subject to ψ = 0 on y = 0 ,ψ → x + y → ∞ ,ψ + U A ( y ) = 0 on C , ( u, v ) · ˆ t = u t on C , (5.7)where C is the vortex boundary and ˆ t is the tangent vector on C . Here u t is the tangentialvelocity inside the vortex which is unknown at this stage. The third boundary conditionis the no-normal flow condition which states that the vortex boundary is a streamline ofthe total streamfunction Ψ = ψ + U A ( y ). Note that this solution is only valid outside thevortex since there are no streamlines which leave the vortex. Therefore inside the vortexwe must instead impose F . Inside the vortex, ψ satisfies1 H ∂ ψ∂x + ∂∂y (cid:20) H ∂ψ∂y (cid:21) + (cid:15) = HF ( ψ + U A ( y )) , (5.8)subject to ψ = 0 on y = 0 ,ψ + U A ( y ) = 0 on C , ( u, v ) · ˆ t = u t on C . (5.9)To obtain a full solution we need to determine the vortex boundary, C and the function, F . This can be achieved by fixing F and then determining C from the requirement that φ + U A = 0 on C and u t is the continuous across C . The typical analytical approachwould be instead to impose a boundary and then use u t from the exterior solution toset a parameter in F (Meleshko & van Heijst 1994; Moffatt 1969), however due to thecomplicated functional dependence on H and nonlinear nature of this problem, there isunlikely to exist a simple expression for the boundary.To proceed we take F ( z ) = − K z, (5.10)inside the vortex. In the case of (cid:15) = 0 and H = const., the system can be solved exactly toobtain the circular Lamb-Chaplygin dipole (Meleshko & van Heijst 1994) vortex solutioncentered at a point on the wall. For the Lamb-Chaplygin dipole we have K = j /a where j ≈ .
83 is the first non-zero root of the Bessel function, J ( z ), and a is the vortexradius. For the case of arbitrary (cid:15) and H we expect K = K ( (cid:15), (cid:104) H (cid:105) , a ) where (cid:104) H (cid:105) denotessome list of parameters of H and does not depend on y . Therefore imposing a value for K sets the size of the vortex as a function of ( (cid:15), (cid:104) H (cid:105) ). We note that the vortex boundary, C , will not necessarily be a semi-circle as in the Lamb-Chaplygin dipole case hence a here ropagation and decay of a shelf vortex (cid:15) (cid:46) H ( y ) varies slowly inside the vortex.We can now seek numerical solutions by choosing a value for K and solving Eqs. (5.6)and (5.8) to obtain the streamfunction, ψ , and hence the vortex size and velocity fields.5.1. Numerical solution
Noting that sgn[ ψ + U A ] = − sgn[ U ] , (5.11)inside the vortex and sgn[ ψ + U A ] = sgn[ U ] , (5.12)outside we can combine Eqs. (5.6) and (5.8) as ∇ ψ − H y H ∂ψ∂y + (cid:15)H = (cid:15)H H ( A − ( ψ/U + A )) θ ( ψ/U + A ) − K H ( ψ + U A ) θ ( − ψ/U − A ) , (5.13)where θ is the heaviside function. Substituting ψ = √ Hφ gives ∇ φ + 12 (cid:34) H yy H − H y H (cid:35) φ = − (cid:15) √ H + (cid:15)H / H (cid:16) A − (cid:16) √ Hφ/U + A (cid:17)(cid:17) θ ( √ Hφ/U + A ) − K H / (cid:16) √ Hφ + U A (cid:17) θ ( −√ Hφ/U − A ) , (5.14)and noting that (cid:15)H / H (cid:16) A − (cid:16) √ Hφ/U + A (cid:17)(cid:17) = (cid:15) √ H − (cid:15)H y φU H + O ( φ ) , (5.15)we may split Eq. (5.14) into linear and nonlinear parts as (cid:34) ∇ + (cid:32) K H θ ( a − r ) + (cid:15)H y U H θ ( r − a ) + 12 (cid:34) H yy H − H y H (cid:35)(cid:33)(cid:35) φ = − √ H (cid:2) (cid:15) + K HU A (cid:3) θ ( a − r ) + N ( φ ) , (5.16)where N ( φ ) = − K [ θ ( − ψ/U − A ) − θ ( a − r )] (cid:104) H φ + H / U A (cid:105) + θ ( r − a ) (cid:20) − (cid:15) √ H + (cid:15)H y φU H (cid:21) + θ ( √ Hφ/U + A ) (cid:15)H / H (cid:16) A − (cid:16) √ Hφ/U + A (cid:17)(cid:17) . (5.17)We note that Eq. (5.16) is an exact rearrangement of Eq. (5.14) for all values of a .However, picking a to be close to the average vortex radius will minimise the nonlinearterm, N , and make finding solutions easier numerically. Defining the left-hand side linearoperator as L = ∇ + (cid:32) K H θ ( a − r ) + (cid:15)H y U H θ ( r − a ) + 12 (cid:34) H yy H − H y H (cid:35)(cid:33) , (5.18)and the φ independent term as C = −√ H (cid:2) (cid:15) + K HU A (cid:3) θ ( a − r ) , (5.19)8 M. N. Crowe & E. R. Johnson we may write Eq. (5.16) as L φ = C + N ( φ ) . (5.20)We can now solve Eq. (5.20) using an iterative method. We begin by finding the linearsolution, φ , satisfying L φ = C, (5.21)by numerically inverting L and imposing boundary conditions of φ = 0 on the boundariesof the numerical domain. Using this linear solution as our initial guess we may now take φ n = L − [ C + N ( φ n − )] , (5.22)and iterate until the domain averaged difference between consecutive φ n is small, (cid:90) D | φ n − φ n − | d A < δ, (5.23)for some δ . Picking a value of a close to the size of the vortex minimises the nonlinearterm N and leads to faster and more consistent convergence. It is sufficient to pick initialradius a using the Lamb-Chaplygin dipole value of a = j /K for small β and (cid:15) . For largerparameter values we can use gradually increase (cid:15) or β while adjusting a to match theobserved vortex size from the previous parameter values. if convergence from the linearsolution is slow or fails, we can use a parameter continuation approach by taking thenonlinear solution for slightly smaller (cid:15) or β as our initial guess, φ .5.2. Results
To illustrate our results we consider H given by Eq. (2.7). The cross-sectional area, A ,is now given by A ( y ) = (cid:40) β [exp[ βy ] − , y (cid:54) D, β [exp[ βD ] −
1] + ( y − D ) exp[ βD ] , y (cid:62) D, (5.24)which can be easily inverted for A − ( z ).We will consider here the case where the vortex radius is less than the shelf width, a < D , however our theory and numerical method is also valid for a > D . Therefore thelinear operator, L , can be written as L = ∇ + K exp[2 βy ] θ ( a − r ) + (cid:18) (cid:15)βU θ ( r − a ) − β (cid:19) θ ( D − y ) , (5.25)which limits to the Linear operator for the Lamb-Chaplygin dipole problem L = ∇ + K θ ( a − r ) , (5.26)in the case of small (cid:15), β . The terms C and N ( φ ) from Eq. (5.20) can be similarly expressedin terms of β and D .Figs. 8 and 9 shows our nonlinear solutions for ψ for a range of parameters. Fig. 8 showsvortices which travel in the opposite direction to the wave field ( U = −
1) for parametervalues ( (cid:15), β ) = { . , , } × { . , } while Fig. 9 shows vortices which travel faster thanthe fastest topographic mode for (cid:15) = 0 . β ∈ { . , } . The solutions are calculatedon the numerical domain ( x, y ) ∈ [ − . , . × [0 , .
2] using (
N x, N y ) = (2048 , δ = 10 − in Eq. (5.23). Weplot the streamfunction in the frame of the vortex, ψ + U A ( y ), so the streamline of height0 denotes the vortex boundary. While our system depends on 5 parameters, (cid:15) , β , D , U and K , we can set | U | = 1 and fix K . Setting | U | = 1 is equivalent to setting the velocity ropagation and decay of a shelf vortex -2 -1 0 1 2 x y -2-1.5-1-0.500.5 (a) -2 -1 0 1 2 x y -6-4-20 (b) -2 -1 0 1 2 x y -2-1.5-1-0.500.5 (c) -2 -1 0 1 2 x y -6-4-20 (d) -2 -1 0 1 2 x y -2-1.5-1-0.50 (e) -2 -1 0 1 2 x y -6-4-20 (f) Figure 8: Plots of the streamfunction in the vortex frame, ψ + U A ( y ), for 6 pairs ofparameters, ( (cid:15), β ), given by (a) (0 . , . . , , . , , . , K = j , U = − D = 12 .
5. Dependence on D is weak and hence not shown here. The dotted line shows the vortex boundary for theLamb-Chaplygin dipole case of ( (cid:15), β ) = (0 ,
0) using the same value of K .scale used for calculating the inverse Rossby number, (cid:15) , whereas fixing K determines thesize of the vortex and hence we can measure the slope, β , and shelf width, D , in unitsof the Lamb-Chaplygin dipole radius, a = j /K . We therefore vary only (cid:15) , β and D andnote that the dependence of the vortex structure on D is weak and not shown here. Wenote, however, that D can play an important role in setting the speed of the modes andhence is chosen such that the vortices in Fig. 9 do not generate topographic modes. Forthe decaying vortex problem where energy is emitted towards the edge of the shelf weexpect that D will play a more important role.From Figs. 8 and 9 we observe that the effect of increasing β is to reduce the vortexsize such that the vortex remains approximately the same shape. This corresponds torequiring a smaller value of the internal wavenumber, K , to get a vortex of approximateradius equal to 1, and hence corresponds to a reduction of maximum vorticity whencompared to an Lamb-Chaplygin dipole of the same size. Increasing (cid:15) has the effectof squashing the vortex against the wall. This is to be expected as the large (cid:15) limitcorresponds to quasi-geostrophic (QG) dynamics in which the flow must follow contours0 M. N. Crowe & E. R. Johnson -2 -1 0 1 2 x y -0.500.511.52 (a) -2 -1 0 1 2 x y (b) Figure 9: Plots of the streamfunction in the vortex frame, ψ + U A ( y ), for (cid:15) = 0 . U = 1 and β = 0 . β = 1 (b). Solutions are shown for K = j and D = 12 . (cid:15), β ) = (0 ,
0) using the same value of K . These solution correspond to vortices travellingfaster than the fastest topographic mode. While the dependence of the vortex structureon D is weak, the value of D = 12 . β , the solution must adjust such that the rotating fluid inside the vortex anddeflected fluid outside the vortex move predominantly in the x direction.
6. Discussion and conclusions
Here we have considered the evolution of a vortex moving along a shelf. To allow forthe calculation of analytical results in the limit of shallow slope and slow rotation, wetake our vortex to be contained within an approximately semi-circular region against thewall. Vortices of this form therefore limit to one half of the Lamb-Chaplygin dipole withthe other half corresponding to the vortex image.Since shelf system admit shelf wave solutions, we began by determining the speedand structure of these modes. For positive rotation rate, these shelf waves move withthe coastal boundary on the right as expected for coastal trapped waves. The finiteshelf width acts to discretise the modes leading to a countable set of wave solutions.The alongshore and offshore wavenumbers can be determined numerically by solving atranscendental equation and the frequency and phase speed of each mode can then bedetermined.If the speed of the vortex matches the phase speed of any shelf wave modes, we expectthe vortex to excite these modes generating a wave wake. Using a Fourier transformapproach, we have determined the far-field amplitude of the wave wake and providedanalytic predictions for the number of modes generated as a function of the vortex speed,shelf parameters and rotation rate. We observe that a slower vortex will match the phasespeeds of a greater number of waves and will hence generate a higher number of differentmodes. The group velocity of each modes is less than or equal to its phase velocity henceall modes will be emitted behind the vortex and we expect no upstream wave signature.The generation of waves corresponds to a flux of energy from the vortex to the wavefield with this flux resulting in the slow decay of the vortex. We have determined theleading order wave energy flux using our far field wave solution and by equating this fluxto the loss of vortex energy we can describe the vortex decay. This decay is shown to beproportional square of the vortex dipole strength and to have a complicated dependenceon the vortex speed, rotation rate and shelf slope. The vortex slows as it loses energy, ropagation and decay of a shelf vortex
Supplementary material.
A supplementary movie is available.
Funding.
This work was funded by the UK Natural Environment Research Councilunder grant number NE/S009922/1.
Declaration of Interests.
The authors report no conflict of interest.
REFERENCESBretherton, F. P.
J. Fluid Mech. (3), 545–570. Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D. & Brown, B. P.
Phys. Rev. Res. ,023068. Crowe, M. N. & Johnson, E. R.
J. Fluid Mech. , A22. M. N. Crowe & E. R. Johnson
Crowe, M. N., Kemp, C. J. D. & Johnson, E. R.
J. Fluid Mech. (in review) . Deremble, B., Johnson, E.R. & Dewar, W.K.
Ocean Modelling , 1 – 12.
Dewar, W. K., Berloff, P. & Hogg, A. McC.
J. Mar. Res. (4-5), 501–522. Dewar, W. K. & Hogg, A. McC.
Ocean Modelling (1), 1 – 13. Flierl, G. R. & Haines, K
Phys.Fluids (10), 3487–3497. Fraenkel, L. E.
Proc. R. Soc. A. (1195), 506–526.
Hill, M. J. M.
Phil. Trans. R. Soc. (A.) , 213–245.
Hogg, A. McC., Dewar, W. K., Berloff, P. & Ward, M. L.
J. Fluid Mech. ,194–208.
Isern-Fontanet, J., Garc´ıa-Ladona, E. & Font, J. .2006 Vortices of the mediterraneansea: An altimetric perspective.
J. Phys. Oceanogr. (1), 87 – 103. Johnson, E. R.
J. Fluid Mech. , 69–76.
Johnson, E. R. & Crowe, M. N.
J. Fluid Mech. (in review) . Johnson, E. R. & Rodney, J. T.
ContinentalShelf Res. (14), 1481 – 1489. LeBlond, P. H. & Mysak, L. A.
Waves in the Ocean . Elsevier.
Lighthill, M. J.
J. Fluid Mech. (4), 725–752. Long, R. R.
J. Met. , 197–203. Machicoane, N., Labarre, V., Voisin, B., Moisy, F. & Cortet, P.-P.
Phys. Rev. Fluids , 034801. Meleshko, V. V. & van Heijst, G. J. F.
J. Fluid Mech. , 157–182.
Moffatt, H. K.
J. Fluid Mech. (1),117–129. Penduff, T., Juza, M., Barnier, B., Zika, J., Dewar, W. K., Treguier, A.-M., Molines,J.-M. & Audiffren, N.
J. Climate24