- Email: [email protected]
Numerical studies of barotropic modons
Dynamics o f Atmospheres and Oceans, 5 (1981) 219--238
219
Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands
NUMERICAL STUDIES OF BAROTROPIC MODONS
JAMES C. McWILLIAMS 1, GLENN R. F L I E R L 2, VITALY D. LARICHEV 3 and G R E G O R Y M. REZNIK 3
1 National Center for Atmospheric Research, Boulder, CO 80307 (U.S.A.) 2 Department o f Meteorology, Massachusetts Institute of Technology, Cam bridge, MA 02139 (U.S.A.) 3 p.p. Shirsov Institute of Oceanology, Academy o f Sciences of the U.S.S.R., Moscow (U.S.S.R.) (Received July 21, 1980; accepted November 17, 1980)
ABSTRACT McWilliams, J.C., Flierl, G.R., Larichev, V.D. and Reznik, G.M., 1981. Numerical studies of barotropic modons. Dyn. Atmos. Oceans, 5: 219--238. Numerical solutions of barotropic modons are examined to assess the accuracy with which they can be calculated, their behavior under the influence of dissipation, and their resistance to perturbations.
1. INTRODUCTION
Modons * are exact, permanent-form, uniformly propagating solutions to the quasigeostrophic potential vorticity equation, in which there is a balance between linear Rossby wave dispersion and non-linear advective steepening in a manner analogous to soliton balances in other physical circumstances (Zabusky and Kruskal, 1965; Scott et al., 1973). The basic theory of modons has been presented by Stern (1975), Larichev and Reznik (1976), Berestov (1979) and Flierl et al. (1980). This paper reports numerical calculations relating to various properties of barotropic modons. In particular, the calculations illustrate: (i) the accuracy attainable by standard numerical techniques, as a function of the degree of resolution, in simulating the analytically known behavior of inviscid modon propagation and persistence;
* There persists a disagreement among the authors as to the appropriate name for the solutions under consideration (e.g. the Introduction of Flierl et al., 1980). Some favor " m o d o n " , w h i c h can be perceived as indicating r o b u s t n e s s and relevance to geophysical flows, and others the more neutral solitary eddy solution (SES). For this paper, we adopt the former and here simply record our disgreement. 0377-0265/81/0000--0000/$02.50 © 1981 Elsevier Scientific Publishing Company
220
(ii) the influence of various types of m o m e n t u m dissipation on m o d o n evolution; and (iii) the resistance of modons to perturbations with various amplitudes and scale contents (in the presence of weak dissipation). 2. BASIC E Q U A T I O N S
We shall restrict our attention to barotropic flow (i.e. where motions are purely horizontal), in which the governing potential vorticity equation is ~t + J ( ~ , ~) + ~x = F
(1)
where J denotes the two-dimensional Jacobian operator, (12)
~" = V2~//
is the relative vorticity and F = (--I)" K, V 2 ( " - : ) ~ "
(3)
is the rate of vorticity dissipation. F assumes different functional forms for n = I, 2 and 3; these are referred to as linear, Newtonian and biharmonic dissipation laws respectively. Equations 1--3 are written in non-dimensional form, where the horizontal coordinates have been scaled by the m o d o n radius a, the time coordinate b y / 3 - ' a - : (~ is the y derivative of the Coriolis frequency), the stream function ~ by ~a 3 , the velocities by f]a2 and the relative vorticity by f~a. As a consequence the ratio of K, and the corresponding dimensional dissipation coefficient is ~ - 1 a : - 2 " . For F = 0, the m o d o n solution to eq. 1 can be written as
k~c~
1 + ~ -c r,
r
(4)
@(r, 0) = c sin 0!
| K1(r/x/~)
r>
[KI(1/V~) '
1
[J:(kr) | Jl(k)
'
r< 1
~(r, 0) = --sin 0 |~ K l ( r / v ~ )
(5) r>:
where J: and K: are Bessel functions of order one. In (4) and (5), r and 0 are polar coordinates in a frame of reference moving in the x-direction with speed c. These formulas describe a family of solutions, spanned as c varies from zero to positive infinity, with an interior region wavenumber k(c) determined from J2(k) K2(l/xfc) (6) k J l ( k ) - %/c K:(I/x/~)
221
I0
(a)
0.0
~ ~
-1.6i
0.0,
1 . 6
0 I0
(b) ~;
~ _ 1 6 t 0.0
161 0.0
0
i 0
x
i i i
I0
Fig. 1. S t r e a m f u n e t i o n ~ and relative vorticity ~ patterns from eqs. 4 and 5 with c = 1. C o n t o u r ii~tervals (CI) are 0.2 and 2 0 respectively The grid resolution is DS = 0.05 (i e as in expt.~5).
222
We shall restrict our attention to the smallest positive solution of (6) and the particular value c = 1: therefore h = 3.9226. The resulting patterns for ~ and are shown in Fig. 1. Hereinafter we shall only show ~ patterns because they are much more sensitive indicators of spatial detail. 3. NUMERICAL TECHNIQUES
Numerical solutions to (1--3) are obtained in a doubly periodic domain of dimension L. We choose L = 10 (i.e. t 0 m o d o n radii) in order that m o d o n solutions in this domain behave similarly to those in an u n b o u n d e d domain; since the m o d o n amplitude in (4) and (5) is exponentially decaying on a unit scale in r for r > 1, the influence of the neighboring modons imposed by periodicity is of the order of e - : ° , or 1 part in 2 × 104. The differential operations in (1--3) are approximated by centered, second-order-accurate, finite differences ("leapfrog" in t i m e ; t h e Arakawa (1966) formulation for the Jacobian advection term). Values of ~ and ~"are c o m p u t e d on a uniform, two-dimensional grid, with horizontal spacing DS. In each direction there are N + 1 grid points spanning [0, L] ; hence D S = L / N . The time step is denoted by Dt and the extent of the integration by T. To diminish consequences of "time splitting", which under some circumstances is an error associated with leapfrog time stepping (Lilly, 1965), the dependent variables at adjacent time steps are averaged at time intervals of approximately At = 0.2, with a subsequent restart of the leapfrog integration regressed in time by (1/2)Dt. Initial values ~o of the m o d o n vorticity are determined at grid points from (5) with the origin of (r, 0) at x = y = L / 2 . To assure periodicity at t = 0, ~-o at y = 0, L is set to zero. A further degree of smoothness near these boundaries is imposed by applying a three-point smoothing operation ~-smooth = (1/4)~oj_: + (1/2)~-0j + (1/4)~-0j+: i,j
(7)
a t j = 2, 3, N - - 1 and N (i a n d j are x and y indices for the grid points, each of which ranges from 1 to N + 1). Initial values of the streamfunction ~i.~ are determined by numerical inversion of (4). 4. RESOLUTION EXPERIMENTS
A first sequence of numerical experiments was performed to illustrate the degree of accuracy as a function of resolution. The parameters for this sequence are listed in Table I. The values of DS range from 0.42 (expt. 1) to 0.05 (expt. 5). Even at the coarsest resolution (expt. 1) an analog of the continuous m o d o n (4--5) occurs in the numerical solutions. A paired structure in vorticity moves eastward with little change in shape or amplitude and with little wave dispersion (Fig. 2). Its average propagation rate, however, is far too slow compared with the continuous value of c = 1; the numerical rate in this case
223 TABLE I i Resolution experiments. F = 0, and the initial c~nditions contain only the modon values i
Expt. no.
N
DS
Dt(X 10 3 )
T
1 2 3 4 5
24 50 100 150 200
0.417 0.200 0.100 0.067 0.050
12.50 6.25 3.13 1.56 1.56
4.97 4.97 4.98 4.99 4.99
is (xc(T) - - x c ( O ) ) / T +- D S / T = 0.50 + 0.08 where xc(t) is the x-position of the vorti¢ity maximum. Also, a patch of small-scale vorticity is left at the origin as the m o d o n propagates away. This presumably represents the residue of an ~djustment between point sampling of the continuous m o d o n in the initial c~nditions and a numerical analog which can be uniformly propagating and shape-preserving under the finitedifference operations. At the finest resolution (expt. 5) the Lpproximation to continuous m o d o n behavior is much more accurate (Fig. 3). The average propagation rate is close to unity, 0.99 + 0.01, and the resic ual vorticity patch has smaller scale and amplitude than for coarser resolutio as. However, the decrease in the residue amplitude, from around 2.0 in e: ;pt. 1 to 0.5 in expt. 5, is not as great as the formal second-order accurac y of the numerical techniques would predict, at least for the resolutions testeq t. Thus, even at relatively fine resolution, there apparently remains a notic~ ~able distinction between a continuous m o d o n and its numerical analog. Th: s is to be expected for solutions with finite-order discontinuities (Orzsag and Jayne, 1974). A useful summary measure of numeri ~al accuracy as a function of resolution is the average propagation rate. Thi~ is plotted in Fig. 4 for expts. 1--5: the accuracy improves as the resolution acreases, the numerical prediction is within 10% of the theoretical value for 5 or more grid intervals per m o d o n diameter, and there is little practical utility in increasing this n u m b e r beyond 30. For this reason, the remaining experiments in this study were performed at a resolution of 30 (i.e. DS = 0.067). This figure can be compared with its counterpart for linear wave propagation l(Fig. 10 of Thompson, 1961) : the resolution requirements for modons are !more severe by approximately a factor of 2. A practical application can be made ~f the preceding result. Elsewhere it has been argued that some atmospheric ~locking events can usefully be interpreted as m o d o n persistence (McWillian~s, 1980). In t h a t interpretation, the m o d o n radius and speed for a January, 1963 block were 1.3 X 106 m and 14
224
t=O.O
]~~~
-07 15.9 0.0.
0.0
~~07
t =4.97
-159
\o.o J
~oo~/22
Fig. 2. C o n t o u r s o f ~ at t = 0 and 4.97 for e x p t . 1 (DS = 0.42); t h e CI is 2 . 0
225
I0
t=O.O
0.0
0 I =
4.99
Fig. 3. C o n t o u r s o f ~" a t t = 0 a n d 4 . 9 9 f o r exl~t. 5 ( D S = 0 . 0 5 ) ; t h e CI is 2.0.
226 TheoreticolIVolue
1.0
i
i
I 2:0
I 30
1
E3 LU LU W 05 nr
LU >
0
0
[ I0
I 40
NUMBER OF GRID INTERVALS PER MODON DIAMETER
Fig. 4. Average propagation speeds as a function of resolution for expts. 1--5. The error bars are calculated from +DS/T.
m s -1 (the l a t t e r is relative t o a m e a n a d v e c t i o n such t h a t the b l o c k r e m a i n e d g e o s t a t i o n a r y ) . A r e s o l u t i o n o f 5 ° in l a t i t u d e and l o n g i t u d e is n o t u n c o m m o n f o r a global w e a t h e r p r e d i c t i o n calculation. A t such a resolution, t h e r e are a b o u t 5 grid intervals per m o d o n d i a m e t e r : h e n c e , f r o m Fig. 4, an e r r o r in speed o f 50%, or 7 m s -1 . Over the course o f a week, which is n o t long c o m p a r e d with the lifetime o f a block, an actually s t a t i o n a r y p a t t e r n c o u l d be m i s f o r e c a s t as m o v i n g 4 0 0 0 km, the w i d t h o f the N o r t h Atlantic Ocean. 5. DISSIPATION EXPERIMENTS T h e t h r e e d i f f e r e n t f u n c t i o n a l f o r m s f o r dissipation in (3) allow comparisons o f e f f e c t s peculiar t o a p a r t i c u l a r f o r m and those which are plausibly general for a b r o a d class o f f o r m s . This d i s t i n c t i o n is o f s o m e i m p o r t a n c e for geophysical i n t e r p r e t a t i o n s o f m o d o n b e h a v i o r because the true n a t u r e o f geophysical dissipation (i.e. m o m e n t u m loss f r o m large-scale flow structures) is n o t k n o w n . A list o f the n u m e r i c a l p a r a m e t e r s f o r the s e q u e n c e o f dissip a t i o n e x p e r i m e n t s is given in Table II. T h e first t h r e e e x p e r i m e n t s were c h o s e n so t h a t a p p r o x i m a t e l y one e-folding d e c a y in the m a x i m u m v o r t i c i t y value o c c u r s in At = 5 t i m e units for each o f the t h r e e n values. F o r n = 1, linear dissipation, the d e c a y is an a p p r o x i m a t e l y shape-preserving e x p o n e n tial f u n c t i o n o f t i m e (it w o u l d be so e x a c t l y if the t e r m s ~x and J(qJ, ~) were a b s e n t f r o m eq. 1), so t h a t the a p p r o p r i a t e value f o r the c o e f f i c i e n t K1 can be d e t e r m i n e d f r o m e--Kl ~t = e--i which yields K I = 0.2. For n = 2 and 3 the decay law is not so simple, and trial and error was required to select the K, value which yields the desired bulk decay. The accuracy of this selection is indicated in the fourth column of Table II (i.e. about +1%).
227 i
TABLE II
Ji
Dissipation experiments. N = 150, DS = 0.067, Dt = 1.56 x 10-3, and the initial conditions contain only the modon values I Expt. no.
n
Kn
(x I y)
n~ax ~'(x, y, t = 4.98) _ max ~(x,y, t = 0)
T
6 7 8 9
1 2 3 2
0.2 0.014 0.001 0.001
01373 ' 0.1370 0.i368 0.i925
4.98 14.94 4.98 4.98
F o r the m o d o n s o l u t i o n (5)
v2(n-1)~/~ = ( - k 2 ) . - 1 f o r r < 1 (which is the region c o n t a i n i n g Ithe ~" m a x i m u m ) . Thus, u n d e r the a s s u m p t i o n o f shape p r e s e r v a t i o n and e x p o n e n t i a l d e c a y , we m i g h t estimate the dissipation c o e f f i c i e n t s a s / ~ n , w h e r e e x p [ - - k 2 ( ~ - l ) k ~ A t ] = e-1 T h e values o b t a i n e d f r o m this f o r m u l a ate k l = 0.2, K2 = 0 . 0 1 3 a n d / { s = 0 . 0 0 0 8 . C o m p a r i n g these e s t i m a t e d c o e f f i c i e n t s with those in Table II shows t h a t the first is an a c c u r a t e estimate, whlle the o t h e r t w o are u n d e r e s t i m a t e s (i.e. the h i g h e r - o r d e r dissipation processes are less e f f i c i e n t t h a n a shape-preserving a r g u m e n t w o u l d indicate). A partial e x p l a n a t i o n o f this e f f e c t is ~hat, in a d d i t i o n t o the p r e c e d i n g ratio v2(n--1)~'/~
" =
C 1-n
for r > 1. Thus, in the o u t e r region with! the same assumptions as above = ~ ( t = O) e x p [ ( - - 1 ) " g ~ c l - ' t ]
.
Since c < k z f o r o u r s o l u t i o n , dissipation is thus locally less e f f i c i e n t in the o u t e r region t h a n for r < 1 w h e n e v e r n :> 1, and an integral dissipation coefficient m u s t be larger t h a n Kn t o a c c o m )lish a specified rate o f decay. This e f f e c t s h o u l d increase with n and is furt: ~er e n h a n c e d f o r even values o f n, w h e r e the d e c a y rate is negative in the c u t e r region. We e x p e c t t h a t an integral r a t h e r t h a n a p o i n t dissipation rate s m o r e relevant d u e t o fl and n o n linear a d j u s t m e n t s (i.e. r e d i s t r i b u t i o n s ) n the p r e s e n c e o f dissipation in o r d e r t o k e e p t h e f l o w s t r u c t u r e close t o a m( don. Even p u r e diffusion c a n n o t be globally shape-preserving ( f o r n = 1) du~ i t o t h e i n f l u e n c e o f t h e fl and n o n linear terms. An alternative e x p l a n a t i o n o f Kn > K~ for n > i is t h a t h i g h e r - o r d e r dissi-
228
pation acts selectively on the smaller scales and, once their amplitude has been diminished, there is less in the flow configuration upon which the process can act efficiently. This argument supposes that there is n o t an efficient non-linear transfer to replenish the depleted energy in the smaller scales such as occurs in a turbulent cascade; in its inviscid form a m o d o n is non-cascading. Figure 5 shows a comparison of meridional sections of ~ for the three dissipation forms after an integration time of At ~ 5. For each form, the maxim u m value in ~ is nearly identical, though greatly reduced from the initial m a x i m u m (Fig. 1). However, the size of the m o d o n expands with larger n. This is plausible since higher-order dissipation forms are relatively more confined in their damping effects to small scales and, after the action of such dissipation, the remaining structures are relatively richer in their large-scale content. Alternatively, we can characterize F for n > 1 as diffusion, and Fig. 5 shows an increase in size with n as the consequence of diffusive spreading of gradients. The decay of the m a x i m u m value of ~ is shown in Fig. 6, again for the three values of n. The preceding argument is further supported by the fact t h a t the decay rate is more rapid for smaller n. Streamfunction, by being a second integral, is a larger-scale structure than vorticity (Fig. 1), and the more small-scale-selective dissipations (with larger n) are less efficient in making ~ decay when the ~ decay rates have been made comparable. As dissipation causes the m o d o n amplitude to diminish in time, the eastward propagation rate also decreases. This is consistent with the m o d o n dispersion relation: as the vorticity m a x i m u m decreases for a fixed m o d o n
i
'
I
i
I
r
I
I
I
I
I
I
1
I
I
'
I
i
....... n=:3
,
I0 ~
,
l
l
l
l
l
l
l
l
--
- - n = 2
6 4
~
2
~
0
~
~~ oo
-2 -4 -6 -8 I
-I0 0
I
2
L
I
3
i
I
4
i
I
5
I
I
6
i
I
7
i
I
8
t
I
9
)-I
J
I0
I
I
1.0
I
2.0
I
I
3.0
L
I
4.0
L
,5.0
Fig. 5. M e r i d i o n a l profiles o f ~ a t x = x c ( t ) a n d t ~ 5 for expts. 6--8. n is a n i n d e x o f t h e f u n c t i o n a l f o r m o f t h e dissipation. Fig. 6. M a x i m u m values o f ~ as a f u n c t i o n o f t i m e for expts. 6--8.
229
radius, c also decreases (because k(c) inc teases from eq. 6 such that the magnitude of J l ( k ) increases in the denomin ator in eq. 5). Physically this p h e n o m e n o n is connected with a decrea ~e in the swirling velocity in each of the separate m o d o n vortex centers and concomitant decrease in the strength of advective pushing towards t~ e east of each vortex center by the field of the other. The remarkable resuli from the set of expts. 6--8 is that the nature of the Xc decrease is identical for the three values of n. This is shown in Fig. 7a. The effect of dissipatign on m o d o n zonal propagation is independent of the form of the dissipatipn and only depends upon its rate, here chosen to produce identical bulk decay in ~"for each value of n. Figure 7b shows the history of Ayc(tl, the meridional separation of the opposite-signed extrema in ~'. As discussedabove, Ay c increases more rapidly as n increases. The sharp changes in the Aylc curves are an artifact of the discrete analysis. The smallest change A y c tan make for a meridionally antisymmetric solution is 2DS, and Fig. 7 is based upon a time sampling of the numerical solution of At = 0.1; hence, the ~lope of the abrupt changes in Fig. 7b is 2 D S / A t = 1.3. Figures 5--7 only span a time interval~n which xc > 0 and we can interpret the evolution as a m o d o n of declining amplitude and Xc. There are, however, threshold values for these quantities implied by the m o d o n theory; namely max ~= 1.7 (x,y) .,~c = 0
for unit m o d o n radius. For any value of the dissipation coefficient, after
I0
.
(a]r
/ I
i
I
'
I
i
II
2.0
I /
t
1.8
\~ ~s6~
8
-
1.4
6
1.2
i
I
I
I
!
""if
/-
jr-
n=l
0.8
---°=' - -
I
•." z./
1.0
!}
t
1.6
7
4
I
(b)
0.6 ~ 04-
n=2
n=3
O2 0
I 0
1.0
L
I ~' ,0
i
I 3.0
i
I 4.0
i ,50
0
I
0
1.0
i
I
2.0
t
I
3.0
t
I 4.0
1 5.0
Fig. 7. Time variations of (a) the zonal positio] Lof the m o d o n e x t r e m a in ~', Xc(t), and (b) the meridional separation of the extrema, Ay c, t), for expts. 6--8.
230
J
I
I
I /
I
'
I
l
_
o+ (3O +.
~J
i
0
I o
i
I o
i
I O
i
I o
o o
(l) ~*
o
'
L"f'-
I
'
I
'
I
i
I
i
I
l
I
I
I
I
0
231
OO 0.0 -0.8
.0 0.0
,~~~-7.1 0.0
0.8
0.0
cb~ ,-,2
\
~oo~ •
-0.~
O0
\oOo~
00
\I
~' 0.0
o,
...oo-
~o,"
x Fig. 9. C o n t o u r s o f ~ at times in t w o different regimes. Panel (a) is in the m o d o n regime, and panel (b) is in the Rossby wave regime. The CI are 0.8 and 0.2 respectively.
232
sufficient time the threshold will be reached and a m o d o n interpretation must be abandoned. Experiment 7 with n = 2 (Newtonian dissipation) was continued to T ~ 15, well past the threshold. Histories of max ~"and Xc are plotted in Fig. 8. The eastward propagation ceases around t = 8 where max is about 3. This is larger than the above threshold value for unit radius, but we recall from Figs. 5 and 7b t h a t the radius has become larger than one and, from the non-dimensionalization, that m o d o n vorticity scales linearly with radius. For times greater than 8, propagation becomes westward, and the decline in max ~"continues. This behavior is crudely t h a t of dissipating, dispersing Rossby waves. The dispersion at later times is illustrated in Fig. 9. Prior to the threshold (panel (a)), there is little dispersion; subsequently (panel (b)), there is a considerable amount. 6. P E R T U R B A T I O N E X P E R I M E N T S
A sequence of experiments has also been performed to illustrate the resistance of modons to perturbations (Table III). In this sequence, the initial vorticity at grid points is the sum of the m o d o n vorticity (as described previously) and a perturbation eZij. Zij is evaluated, at grid points (xi, Yj), from M u
Z(x, y) = Co
~
a,m cos(2~nx/L + O,m) cos(2~my/L + O,m)
(8)
n, m =M1
In (8), the arrays a,m, Onto and ~,m are chosen from a computer random number generator which approximates a uniform probability distribution with mean zero and a range of either [ - x / 1 2 , +x/12] (for a,m, so that a2m = 1)
T A B L E III P e r t u r b a t i o n e x p e r i m e n t s . N = 150, DS = 0.067, D t = 1.56 x 10 - 3 , n = 2 and K 2 = 0.001 E x p t . no.
Ml
Mu
e
RMS t~ *
T
Modon survival?
10 11 12 13 ** 14 15 16 17
1 1 1 1 1 1 1 6
5 5 5 5 5 25 25 25
0.1 0.4 1.0 1.0 4.0 1.0 4.0 4.0
0.04 0.15 0.37 0.38 1.50 0.08 0.31 0.04
4.98 4.98 14.94 9.96 4.98 9.96 9.96 14.94
Yes Yes No No No Yes No Yes
* RMS ~ is the r o o t m e a n square p e r t u r b a t i o n s t r e a m f u n c t i o n , b a s e d u p o n an area avera g e ; i t is a s o l u t i o n o f V2t~ = eZ, w h e r e Z is given b y eqs. 8 and 9. ** E x p e r i m e n t 13 is identical t o e x p t . 12 e x c e p t t h a t it is based u p o n an alternative set o f r a n d o m n u m b e r s for t h e p e r t u r b a t i o n a m p l i t u d e s and phases.
233
or [--w, ~r] (for 0nm and ¢,m). The constant Co in (8) is chosen such that the root mean square (RMS) perturbation vorticity is unity when e = 1; that is N
Zi2j / N 2 = 1
(9)
i,j = 1
The sequence in Table III represents different perturbation amplitudes (e) and scale contents (MI, Mu). We can identify a modon diameter with a perturbation wavelength L / M for wavenumber M = 5. Thus, there are perturbations with scales as large as or larger than the modon scale (expts. 10--14) or with scales smaller (expt. 17) or with scales spanning the modon scale (expts. 15, 16). The perturbation experiments have a weak dissipation (equivalent to that of expt. 9). The value of K2 is small enough to alter only slightly the modon propagation and persistence (in expt. 9 max ~"and Xc are within 10% of their inviscid values after At ~ 5), yet it provides a sink for enstrophy which can accumulate on the grid scale when e is large enough that a twodimensional turbulent cascade occurs. For weak perturbations (small e) modons are unaffected in their behavior, at least for the integration times tested. For example, the amplitude and central position of the modon in expt. 10 are the same at t = 5 as they are in expt. 9, to within an uncertainty no larger than the perturbation amplitude.
i
i
\ s
\
\
D
D
\ \
~struction~ \ \
¢
S
I
\
\ \ \
D \
~urvival~
Ol
o
\ S
r i
i 2
P 3
5 (M~MI) -I/2
Fig. 10. A regime diagram for the survival (S) or d e s t r u c t i o n ( D ) o f a barotropic m o d o n in the presence o f a r a n d o m perturbation. This diagram is based u p o n the results o f expts. 1 0 - - 1 7 . e is the area average RMS v o r t i c i t y a m p l i t u d e o f the perturbation ( t o be c o m pared w i t h a m o d o n v o r t i c i t y m a x i m u m o f 1 6 . 1 ) , and the perturbation scale is a geometric m e a n scale ( n o r m a l i z e d so that u n i t y is the m o d o n scale).
234 M=, M= - 1,5
rms = 0.37
MI,M= - 6 , 2 5
C . I . • 0.1 M,, M, • I, 25
rm$ , 0 . 0 1
rms • 0.08
C.1. • 0 . 0 3
C.1. • 0.01
o_
Fig. 11. Perturbation streamfunetion ~ ( x , y) for three sets of wavenumber intervals (Mu, M1) obtained by solving V2~ = Z, where Z is determined from eq. 8. For each panel the (area average) RMS value for • is listed.
For sufficiently strong perturbations, however, the m o d o n is destroyed as a separate entity by participation in a turbulent cascade. We have used a subjective criterion for m o d o n survival or destruction, y e t in our solutions the issue does n o t appear to be a subtle one because the modon, while surviving, is readily identified (the first panels in Figs. 12, 13). The critical perturbation amplitude required for destruction is a function
235 tmO
1~ 0.4
t =0.8
y
3(
Fig. 12. Vortieity patterns for an example of rapid modon destruction (expt. 14). The CI is 2.0.
of the scale c o n t e n t of the perturbation. Figure 10 is a regime diagram in p e r t u r b a t i o n vorticity amplitude and geometric mean horizontal scale. The critical amplitude decreases as the scale increases. Stated alternatively, vorticity on larger scales is more efficient in destroying m o d o n s than vorticity on smaller scales. Perturbation s t r e a m f u n c t i o n patterns f or e = 1 and various scales are shown in Fig. 11. The RMS amplitude o f • increases with pert urbat i on scale
236
tzl
t :10
Fig. 13. Vorticity patterns for an example of slow modon destruction (expt. 12}. The CI is 1.0.
237 approximately as its square, as expected from a second integral of vorticity. Thus, the boundary between m o d o n survival and destruction is less scaledependent when the perturbation amplitude is measured in terms of stream function rather than vorticity. This boundary occurs at an RMS ( e ~ ) value around 10--20% of the m o d o n e x t r e m u m in streamfunction. The mechanism for m o d o n destruction by perturbation appears to be advective shearing o u t of the m o d o n vorticity contours (which are material lines approximately) b y the large scale perturbations. Identifying the latter with the "energetic eddies", this mechanism is known to be active in twodimensional turbulence (Kraichnan, 1971). This process is indicated b y Fig. 12, which is an example of rapid m o d o n destruction. For smaller values of e, slightly above the critical value, the destruction rate is lower, and the shearing o u t of m o d o n vorticity contours seems to be selectively occurring at the edges of the vorticity centers and subsequently eroding inward. An early consequence of this is a separation of the two m o d o n centers, so that they cease to interact strongly and can begin to move westward and meridionally and otherwise act as isolated vortices (McWilliams and Flierl, 1979). This is illustrated in Fig. 13, which is an example of slow m o d o n destruction. The behavior shown in Fig. 13, which is based on expt. 12, is qualitatively the same in expts. 13 and 16 as well. Thus, the nature of the m o d o n destruction process in the slightly supercritical regime seems moderately independent of either scale content or the particular rand o m realization of the perturbation field. 7. CONCLUSIONS Numerical solutions of the barotropic, ~-plane potential vorticity equation have been examined to address the three aspects of m o d o n behavior listed in the introduction. In summary the results are as follows. (i) Modons can be successfully calculated b y standard numerical techniques if the resolution scales in space and time are sufficiently small. In particular, a b o u t 20 grid points per m o d o n diameter are required to obtain greater than 95% accuracy in the bulk propagation rate using second-order finite~lifference techniques. (ii) Under the influence of m o m e n t u m dissipation, modons decrease in amplitude, reduce their zonal propagation rate and expand their meridional scale. The first two processes occur in ways which are insensitive to gross aspects of the nature of the dissipation, and the third is simply related to the order af the dissipation law. After sufficient amplitude decay the m o d o n structures make a transition to a dispersive Rossby wave regime. While still within the m o d o n regime the decline in amplitude and speed crudely follows a m o d o n dispersion curve. (iii) Modon$ are resistant to perturbations of small amplitude and are destroyed b y perturbations of moderate amplitude. The critical amplitude for destruction is dependent u p o n the scale c o n t e n t of the perturbation, in a
238
manner consistent with advective shearing on larger than m o d o n scale being the dominant destructive mechanism.
ACKNOWLEDGMENTS
The numerical calculations were performed at the National Center for Atmospheric Research, which is sponsored by the National Science Foundation. Programming was done b y Dr. Julianna H.S. Chow. The international collaboration which this paper represents arose under the auspices of the POLYMODE Program.
REFERENCES Arakawa, A., 1966. Computational design for long term numerical integration of the equations of fluid motion: Two dimensional incompressible flow. Part I. J. Comput. Phys., 1 : 119--143. Berestov, A.L., 1979. Solitary Rossby waves, Izv. Acad. Sci., U.S.S.R., Atmos. Oceanic Phys., 15 : 648--654. Flierl, G.R., Larichev, V.D., McWilliams J.C. and Reznik, G.M., 1980. The dynamics of baroclinic and barotropic solitary eddies. Dyn. Atmos. Oceans, 5: 1--41. Kraichnan, R.H., 1971: Inertial-range transfer in two- and three-dimensional turbulence. J. Fluid Mech., 47: 525--535. Larichev, V. and Reznik, G., 1976. Two-dimensional Rossby soliton: an exact solution. Rep. U.S.S.R., Acad. Sci., 231 (5) (also POLYMODE News, 19). Lilly, D.K., 1965. On the computational stability of numerical solutions of time-dependent non-linear geophysical fluid dynamics problems. Mon. Weather Rev., 9 3 : 1 1 - - 2 6 . McWilliams, J.C., 1980. An application of equivalent modons to atmospheric blocking. Dyn. Atmos. Oceans, 5 : 43--66. McWilliams, J.C. and Flierl, G.R., 1979. On the evolution of isolated, non-linear vortices. J. Phys. Oceanogr., 9: 1155--1182. Orzsag, S.A. and Jayne L.W., 1974. Local errors of approximate solutions to hyperbolic problems. J. Comput. Phys., 14: 93--103. Scott, A.C., Chu, F.Y.F. and McLaughlin, D.W., 1973. The soliton : a new concept in applied science. Proc. IEEE, 61: 1443--1483. Stern, M.E., 1975. Minimal properties of planetary eddies. J. Mar. Res., 33: 1--13. Thompson, P.D., 1961. Numerical Weather Analysis and Prediction. Macmillan, New York, NY, 170 pp. Zabusky, N.J. and Kruskal, M.D., 1965. Interaction of "solitons" in a collisonless plasma and the occurrence of initial states. Phys. Rev. Lett., 15: 240--243.
219
Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands
NUMERICAL STUDIES OF BAROTROPIC MODONS
JAMES C. McWILLIAMS 1, GLENN R. F L I E R L 2, VITALY D. LARICHEV 3 and G R E G O R Y M. REZNIK 3
1 National Center for Atmospheric Research, Boulder, CO 80307 (U.S.A.) 2 Department o f Meteorology, Massachusetts Institute of Technology, Cam bridge, MA 02139 (U.S.A.) 3 p.p. Shirsov Institute of Oceanology, Academy o f Sciences of the U.S.S.R., Moscow (U.S.S.R.) (Received July 21, 1980; accepted November 17, 1980)
ABSTRACT McWilliams, J.C., Flierl, G.R., Larichev, V.D. and Reznik, G.M., 1981. Numerical studies of barotropic modons. Dyn. Atmos. Oceans, 5: 219--238. Numerical solutions of barotropic modons are examined to assess the accuracy with which they can be calculated, their behavior under the influence of dissipation, and their resistance to perturbations.
1. INTRODUCTION
Modons * are exact, permanent-form, uniformly propagating solutions to the quasigeostrophic potential vorticity equation, in which there is a balance between linear Rossby wave dispersion and non-linear advective steepening in a manner analogous to soliton balances in other physical circumstances (Zabusky and Kruskal, 1965; Scott et al., 1973). The basic theory of modons has been presented by Stern (1975), Larichev and Reznik (1976), Berestov (1979) and Flierl et al. (1980). This paper reports numerical calculations relating to various properties of barotropic modons. In particular, the calculations illustrate: (i) the accuracy attainable by standard numerical techniques, as a function of the degree of resolution, in simulating the analytically known behavior of inviscid modon propagation and persistence;
* There persists a disagreement among the authors as to the appropriate name for the solutions under consideration (e.g. the Introduction of Flierl et al., 1980). Some favor " m o d o n " , w h i c h can be perceived as indicating r o b u s t n e s s and relevance to geophysical flows, and others the more neutral solitary eddy solution (SES). For this paper, we adopt the former and here simply record our disgreement. 0377-0265/81/0000--0000/$02.50 © 1981 Elsevier Scientific Publishing Company
220
(ii) the influence of various types of m o m e n t u m dissipation on m o d o n evolution; and (iii) the resistance of modons to perturbations with various amplitudes and scale contents (in the presence of weak dissipation). 2. BASIC E Q U A T I O N S
We shall restrict our attention to barotropic flow (i.e. where motions are purely horizontal), in which the governing potential vorticity equation is ~t + J ( ~ , ~) + ~x = F
(1)
where J denotes the two-dimensional Jacobian operator, (12)
~" = V2~//
is the relative vorticity and F = (--I)" K, V 2 ( " - : ) ~ "
(3)
is the rate of vorticity dissipation. F assumes different functional forms for n = I, 2 and 3; these are referred to as linear, Newtonian and biharmonic dissipation laws respectively. Equations 1--3 are written in non-dimensional form, where the horizontal coordinates have been scaled by the m o d o n radius a, the time coordinate b y / 3 - ' a - : (~ is the y derivative of the Coriolis frequency), the stream function ~ by ~a 3 , the velocities by f]a2 and the relative vorticity by f~a. As a consequence the ratio of K, and the corresponding dimensional dissipation coefficient is ~ - 1 a : - 2 " . For F = 0, the m o d o n solution to eq. 1 can be written as
k~c~
1 + ~ -c r,
r
(4)
@(r, 0) = c sin 0!
| K1(r/x/~)
r>
[KI(1/V~) '
1
[J:(kr) | Jl(k)
'
r< 1
~(r, 0) = --sin 0 |~ K l ( r / v ~ )
(5) r>:
where J: and K: are Bessel functions of order one. In (4) and (5), r and 0 are polar coordinates in a frame of reference moving in the x-direction with speed c. These formulas describe a family of solutions, spanned as c varies from zero to positive infinity, with an interior region wavenumber k(c) determined from J2(k) K2(l/xfc) (6) k J l ( k ) - %/c K:(I/x/~)
221
I0
(a)
0.0
~ ~
-1.6i
0.0,
1 . 6
0 I0
(b) ~;
~ _ 1 6 t 0.0
161 0.0
0
i 0
x
i i i
I0
Fig. 1. S t r e a m f u n e t i o n ~ and relative vorticity ~ patterns from eqs. 4 and 5 with c = 1. C o n t o u r ii~tervals (CI) are 0.2 and 2 0 respectively The grid resolution is DS = 0.05 (i e as in expt.~5).
222
We shall restrict our attention to the smallest positive solution of (6) and the particular value c = 1: therefore h = 3.9226. The resulting patterns for ~ and are shown in Fig. 1. Hereinafter we shall only show ~ patterns because they are much more sensitive indicators of spatial detail. 3. NUMERICAL TECHNIQUES
Numerical solutions to (1--3) are obtained in a doubly periodic domain of dimension L. We choose L = 10 (i.e. t 0 m o d o n radii) in order that m o d o n solutions in this domain behave similarly to those in an u n b o u n d e d domain; since the m o d o n amplitude in (4) and (5) is exponentially decaying on a unit scale in r for r > 1, the influence of the neighboring modons imposed by periodicity is of the order of e - : ° , or 1 part in 2 × 104. The differential operations in (1--3) are approximated by centered, second-order-accurate, finite differences ("leapfrog" in t i m e ; t h e Arakawa (1966) formulation for the Jacobian advection term). Values of ~ and ~"are c o m p u t e d on a uniform, two-dimensional grid, with horizontal spacing DS. In each direction there are N + 1 grid points spanning [0, L] ; hence D S = L / N . The time step is denoted by Dt and the extent of the integration by T. To diminish consequences of "time splitting", which under some circumstances is an error associated with leapfrog time stepping (Lilly, 1965), the dependent variables at adjacent time steps are averaged at time intervals of approximately At = 0.2, with a subsequent restart of the leapfrog integration regressed in time by (1/2)Dt. Initial values ~o of the m o d o n vorticity are determined at grid points from (5) with the origin of (r, 0) at x = y = L / 2 . To assure periodicity at t = 0, ~-o at y = 0, L is set to zero. A further degree of smoothness near these boundaries is imposed by applying a three-point smoothing operation ~-smooth = (1/4)~oj_: + (1/2)~-0j + (1/4)~-0j+: i,j
(7)
a t j = 2, 3, N - - 1 and N (i a n d j are x and y indices for the grid points, each of which ranges from 1 to N + 1). Initial values of the streamfunction ~i.~ are determined by numerical inversion of (4). 4. RESOLUTION EXPERIMENTS
A first sequence of numerical experiments was performed to illustrate the degree of accuracy as a function of resolution. The parameters for this sequence are listed in Table I. The values of DS range from 0.42 (expt. 1) to 0.05 (expt. 5). Even at the coarsest resolution (expt. 1) an analog of the continuous m o d o n (4--5) occurs in the numerical solutions. A paired structure in vorticity moves eastward with little change in shape or amplitude and with little wave dispersion (Fig. 2). Its average propagation rate, however, is far too slow compared with the continuous value of c = 1; the numerical rate in this case
223 TABLE I i Resolution experiments. F = 0, and the initial c~nditions contain only the modon values i
Expt. no.
N
DS
Dt(X 10 3 )
T
1 2 3 4 5
24 50 100 150 200
0.417 0.200 0.100 0.067 0.050
12.50 6.25 3.13 1.56 1.56
4.97 4.97 4.98 4.99 4.99
is (xc(T) - - x c ( O ) ) / T +- D S / T = 0.50 + 0.08 where xc(t) is the x-position of the vorti¢ity maximum. Also, a patch of small-scale vorticity is left at the origin as the m o d o n propagates away. This presumably represents the residue of an ~djustment between point sampling of the continuous m o d o n in the initial c~nditions and a numerical analog which can be uniformly propagating and shape-preserving under the finitedifference operations. At the finest resolution (expt. 5) the Lpproximation to continuous m o d o n behavior is much more accurate (Fig. 3). The average propagation rate is close to unity, 0.99 + 0.01, and the resic ual vorticity patch has smaller scale and amplitude than for coarser resolutio as. However, the decrease in the residue amplitude, from around 2.0 in e: ;pt. 1 to 0.5 in expt. 5, is not as great as the formal second-order accurac y of the numerical techniques would predict, at least for the resolutions testeq t. Thus, even at relatively fine resolution, there apparently remains a notic~ ~able distinction between a continuous m o d o n and its numerical analog. Th: s is to be expected for solutions with finite-order discontinuities (Orzsag and Jayne, 1974). A useful summary measure of numeri ~al accuracy as a function of resolution is the average propagation rate. Thi~ is plotted in Fig. 4 for expts. 1--5: the accuracy improves as the resolution acreases, the numerical prediction is within 10% of the theoretical value for 5 or more grid intervals per m o d o n diameter, and there is little practical utility in increasing this n u m b e r beyond 30. For this reason, the remaining experiments in this study were performed at a resolution of 30 (i.e. DS = 0.067). This figure can be compared with its counterpart for linear wave propagation l(Fig. 10 of Thompson, 1961) : the resolution requirements for modons are !more severe by approximately a factor of 2. A practical application can be made ~f the preceding result. Elsewhere it has been argued that some atmospheric ~locking events can usefully be interpreted as m o d o n persistence (McWillian~s, 1980). In t h a t interpretation, the m o d o n radius and speed for a January, 1963 block were 1.3 X 106 m and 14
224
t=O.O
]~~~
-07 15.9 0.0.
0.0
~~07
t =4.97
-159
\o.o J
~oo~/22
Fig. 2. C o n t o u r s o f ~ at t = 0 and 4.97 for e x p t . 1 (DS = 0.42); t h e CI is 2 . 0
225
I0
t=O.O
0.0
0 I =
4.99
Fig. 3. C o n t o u r s o f ~" a t t = 0 a n d 4 . 9 9 f o r exl~t. 5 ( D S = 0 . 0 5 ) ; t h e CI is 2.0.
226 TheoreticolIVolue
1.0
i
i
I 2:0
I 30
1
E3 LU LU W 05 nr
LU >
0
0
[ I0
I 40
NUMBER OF GRID INTERVALS PER MODON DIAMETER
Fig. 4. Average propagation speeds as a function of resolution for expts. 1--5. The error bars are calculated from +DS/T.
m s -1 (the l a t t e r is relative t o a m e a n a d v e c t i o n such t h a t the b l o c k r e m a i n e d g e o s t a t i o n a r y ) . A r e s o l u t i o n o f 5 ° in l a t i t u d e and l o n g i t u d e is n o t u n c o m m o n f o r a global w e a t h e r p r e d i c t i o n calculation. A t such a resolution, t h e r e are a b o u t 5 grid intervals per m o d o n d i a m e t e r : h e n c e , f r o m Fig. 4, an e r r o r in speed o f 50%, or 7 m s -1 . Over the course o f a week, which is n o t long c o m p a r e d with the lifetime o f a block, an actually s t a t i o n a r y p a t t e r n c o u l d be m i s f o r e c a s t as m o v i n g 4 0 0 0 km, the w i d t h o f the N o r t h Atlantic Ocean. 5. DISSIPATION EXPERIMENTS T h e t h r e e d i f f e r e n t f u n c t i o n a l f o r m s f o r dissipation in (3) allow comparisons o f e f f e c t s peculiar t o a p a r t i c u l a r f o r m and those which are plausibly general for a b r o a d class o f f o r m s . This d i s t i n c t i o n is o f s o m e i m p o r t a n c e for geophysical i n t e r p r e t a t i o n s o f m o d o n b e h a v i o r because the true n a t u r e o f geophysical dissipation (i.e. m o m e n t u m loss f r o m large-scale flow structures) is n o t k n o w n . A list o f the n u m e r i c a l p a r a m e t e r s f o r the s e q u e n c e o f dissip a t i o n e x p e r i m e n t s is given in Table II. T h e first t h r e e e x p e r i m e n t s were c h o s e n so t h a t a p p r o x i m a t e l y one e-folding d e c a y in the m a x i m u m v o r t i c i t y value o c c u r s in At = 5 t i m e units for each o f the t h r e e n values. F o r n = 1, linear dissipation, the d e c a y is an a p p r o x i m a t e l y shape-preserving e x p o n e n tial f u n c t i o n o f t i m e (it w o u l d be so e x a c t l y if the t e r m s ~x and J(qJ, ~) were a b s e n t f r o m eq. 1), so t h a t the a p p r o p r i a t e value f o r the c o e f f i c i e n t K1 can be d e t e r m i n e d f r o m e--Kl ~t = e--i which yields K I = 0.2. For n = 2 and 3 the decay law is not so simple, and trial and error was required to select the K, value which yields the desired bulk decay. The accuracy of this selection is indicated in the fourth column of Table II (i.e. about +1%).
227 i
TABLE II
Ji
Dissipation experiments. N = 150, DS = 0.067, Dt = 1.56 x 10-3, and the initial conditions contain only the modon values I Expt. no.
n
Kn
(x I y)
n~ax ~'(x, y, t = 4.98) _ max ~(x,y, t = 0)
T
6 7 8 9
1 2 3 2
0.2 0.014 0.001 0.001
01373 ' 0.1370 0.i368 0.i925
4.98 14.94 4.98 4.98
F o r the m o d o n s o l u t i o n (5)
v2(n-1)~/~ = ( - k 2 ) . - 1 f o r r < 1 (which is the region c o n t a i n i n g Ithe ~" m a x i m u m ) . Thus, u n d e r the a s s u m p t i o n o f shape p r e s e r v a t i o n and e x p o n e n t i a l d e c a y , we m i g h t estimate the dissipation c o e f f i c i e n t s a s / ~ n , w h e r e e x p [ - - k 2 ( ~ - l ) k ~ A t ] = e-1 T h e values o b t a i n e d f r o m this f o r m u l a ate k l = 0.2, K2 = 0 . 0 1 3 a n d / { s = 0 . 0 0 0 8 . C o m p a r i n g these e s t i m a t e d c o e f f i c i e n t s with those in Table II shows t h a t the first is an a c c u r a t e estimate, whlle the o t h e r t w o are u n d e r e s t i m a t e s (i.e. the h i g h e r - o r d e r dissipation processes are less e f f i c i e n t t h a n a shape-preserving a r g u m e n t w o u l d indicate). A partial e x p l a n a t i o n o f this e f f e c t is ~hat, in a d d i t i o n t o the p r e c e d i n g ratio v2(n--1)~'/~
" =
C 1-n
for r > 1. Thus, in the o u t e r region with! the same assumptions as above = ~ ( t = O) e x p [ ( - - 1 ) " g ~ c l - ' t ]
.
Since c < k z f o r o u r s o l u t i o n , dissipation is thus locally less e f f i c i e n t in the o u t e r region t h a n for r < 1 w h e n e v e r n :> 1, and an integral dissipation coefficient m u s t be larger t h a n Kn t o a c c o m )lish a specified rate o f decay. This e f f e c t s h o u l d increase with n and is furt: ~er e n h a n c e d f o r even values o f n, w h e r e the d e c a y rate is negative in the c u t e r region. We e x p e c t t h a t an integral r a t h e r t h a n a p o i n t dissipation rate s m o r e relevant d u e t o fl and n o n linear a d j u s t m e n t s (i.e. r e d i s t r i b u t i o n s ) n the p r e s e n c e o f dissipation in o r d e r t o k e e p t h e f l o w s t r u c t u r e close t o a m( don. Even p u r e diffusion c a n n o t be globally shape-preserving ( f o r n = 1) du~ i t o t h e i n f l u e n c e o f t h e fl and n o n linear terms. An alternative e x p l a n a t i o n o f Kn > K~ for n > i is t h a t h i g h e r - o r d e r dissi-
228
pation acts selectively on the smaller scales and, once their amplitude has been diminished, there is less in the flow configuration upon which the process can act efficiently. This argument supposes that there is n o t an efficient non-linear transfer to replenish the depleted energy in the smaller scales such as occurs in a turbulent cascade; in its inviscid form a m o d o n is non-cascading. Figure 5 shows a comparison of meridional sections of ~ for the three dissipation forms after an integration time of At ~ 5. For each form, the maxim u m value in ~ is nearly identical, though greatly reduced from the initial m a x i m u m (Fig. 1). However, the size of the m o d o n expands with larger n. This is plausible since higher-order dissipation forms are relatively more confined in their damping effects to small scales and, after the action of such dissipation, the remaining structures are relatively richer in their large-scale content. Alternatively, we can characterize F for n > 1 as diffusion, and Fig. 5 shows an increase in size with n as the consequence of diffusive spreading of gradients. The decay of the m a x i m u m value of ~ is shown in Fig. 6, again for the three values of n. The preceding argument is further supported by the fact t h a t the decay rate is more rapid for smaller n. Streamfunction, by being a second integral, is a larger-scale structure than vorticity (Fig. 1), and the more small-scale-selective dissipations (with larger n) are less efficient in making ~ decay when the ~ decay rates have been made comparable. As dissipation causes the m o d o n amplitude to diminish in time, the eastward propagation rate also decreases. This is consistent with the m o d o n dispersion relation: as the vorticity m a x i m u m decreases for a fixed m o d o n
i
'
I
i
I
r
I
I
I
I
I
I
1
I
I
'
I
i
....... n=:3
,
I0 ~
,
l
l
l
l
l
l
l
l
--
- - n = 2
6 4
~
2
~
0
~
~~ oo
-2 -4 -6 -8 I
-I0 0
I
2
L
I
3
i
I
4
i
I
5
I
I
6
i
I
7
i
I
8
t
I
9
)-I
J
I0
I
I
1.0
I
2.0
I
I
3.0
L
I
4.0
L
,5.0
Fig. 5. M e r i d i o n a l profiles o f ~ a t x = x c ( t ) a n d t ~ 5 for expts. 6--8. n is a n i n d e x o f t h e f u n c t i o n a l f o r m o f t h e dissipation. Fig. 6. M a x i m u m values o f ~ as a f u n c t i o n o f t i m e for expts. 6--8.
229
radius, c also decreases (because k(c) inc teases from eq. 6 such that the magnitude of J l ( k ) increases in the denomin ator in eq. 5). Physically this p h e n o m e n o n is connected with a decrea ~e in the swirling velocity in each of the separate m o d o n vortex centers and concomitant decrease in the strength of advective pushing towards t~ e east of each vortex center by the field of the other. The remarkable resuli from the set of expts. 6--8 is that the nature of the Xc decrease is identical for the three values of n. This is shown in Fig. 7a. The effect of dissipatign on m o d o n zonal propagation is independent of the form of the dissipatipn and only depends upon its rate, here chosen to produce identical bulk decay in ~"for each value of n. Figure 7b shows the history of Ayc(tl, the meridional separation of the opposite-signed extrema in ~'. As discussedabove, Ay c increases more rapidly as n increases. The sharp changes in the Aylc curves are an artifact of the discrete analysis. The smallest change A y c tan make for a meridionally antisymmetric solution is 2DS, and Fig. 7 is based upon a time sampling of the numerical solution of At = 0.1; hence, the ~lope of the abrupt changes in Fig. 7b is 2 D S / A t = 1.3. Figures 5--7 only span a time interval~n which xc > 0 and we can interpret the evolution as a m o d o n of declining amplitude and Xc. There are, however, threshold values for these quantities implied by the m o d o n theory; namely max ~= 1.7 (x,y) .,~c = 0
for unit m o d o n radius. For any value of the dissipation coefficient, after
I0
.
(a]r
/ I
i
I
'
I
i
II
2.0
I /
t
1.8
\~ ~s6~
8
-
1.4
6
1.2
i
I
I
I
!
""if
/-
jr-
n=l
0.8
---°=' - -
I
•." z./
1.0
!}
t
1.6
7
4
I
(b)
0.6 ~ 04-
n=2
n=3
O2 0
I 0
1.0
L
I ~' ,0
i
I 3.0
i
I 4.0
i ,50
0
I
0
1.0
i
I
2.0
t
I
3.0
t
I 4.0
1 5.0
Fig. 7. Time variations of (a) the zonal positio] Lof the m o d o n e x t r e m a in ~', Xc(t), and (b) the meridional separation of the extrema, Ay c, t), for expts. 6--8.
230
J
I
I
I /
I
'
I
l
_
o+ (3O +.
~J
i
0
I o
i
I o
i
I O
i
I o
o o
(l) ~*
o
'
L"f'-
I
'
I
'
I
i
I
i
I
l
I
I
I
I
0
231
OO 0.0 -0.8
.0 0.0
,~~~-7.1 0.0
0.8
0.0
cb~ ,-,2
\
~oo~ •
-0.~
O0
\oOo~
00
\I
~' 0.0
o,
...oo-
~o,"
x Fig. 9. C o n t o u r s o f ~ at times in t w o different regimes. Panel (a) is in the m o d o n regime, and panel (b) is in the Rossby wave regime. The CI are 0.8 and 0.2 respectively.
232
sufficient time the threshold will be reached and a m o d o n interpretation must be abandoned. Experiment 7 with n = 2 (Newtonian dissipation) was continued to T ~ 15, well past the threshold. Histories of max ~"and Xc are plotted in Fig. 8. The eastward propagation ceases around t = 8 where max is about 3. This is larger than the above threshold value for unit radius, but we recall from Figs. 5 and 7b t h a t the radius has become larger than one and, from the non-dimensionalization, that m o d o n vorticity scales linearly with radius. For times greater than 8, propagation becomes westward, and the decline in max ~"continues. This behavior is crudely t h a t of dissipating, dispersing Rossby waves. The dispersion at later times is illustrated in Fig. 9. Prior to the threshold (panel (a)), there is little dispersion; subsequently (panel (b)), there is a considerable amount. 6. P E R T U R B A T I O N E X P E R I M E N T S
A sequence of experiments has also been performed to illustrate the resistance of modons to perturbations (Table III). In this sequence, the initial vorticity at grid points is the sum of the m o d o n vorticity (as described previously) and a perturbation eZij. Zij is evaluated, at grid points (xi, Yj), from M u
Z(x, y) = Co
~
a,m cos(2~nx/L + O,m) cos(2~my/L + O,m)
(8)
n, m =M1
In (8), the arrays a,m, Onto and ~,m are chosen from a computer random number generator which approximates a uniform probability distribution with mean zero and a range of either [ - x / 1 2 , +x/12] (for a,m, so that a2m = 1)
T A B L E III P e r t u r b a t i o n e x p e r i m e n t s . N = 150, DS = 0.067, D t = 1.56 x 10 - 3 , n = 2 and K 2 = 0.001 E x p t . no.
Ml
Mu
e
RMS t~ *
T
Modon survival?
10 11 12 13 ** 14 15 16 17
1 1 1 1 1 1 1 6
5 5 5 5 5 25 25 25
0.1 0.4 1.0 1.0 4.0 1.0 4.0 4.0
0.04 0.15 0.37 0.38 1.50 0.08 0.31 0.04
4.98 4.98 14.94 9.96 4.98 9.96 9.96 14.94
Yes Yes No No No Yes No Yes
* RMS ~ is the r o o t m e a n square p e r t u r b a t i o n s t r e a m f u n c t i o n , b a s e d u p o n an area avera g e ; i t is a s o l u t i o n o f V2t~ = eZ, w h e r e Z is given b y eqs. 8 and 9. ** E x p e r i m e n t 13 is identical t o e x p t . 12 e x c e p t t h a t it is based u p o n an alternative set o f r a n d o m n u m b e r s for t h e p e r t u r b a t i o n a m p l i t u d e s and phases.
233
or [--w, ~r] (for 0nm and ¢,m). The constant Co in (8) is chosen such that the root mean square (RMS) perturbation vorticity is unity when e = 1; that is N
Zi2j / N 2 = 1
(9)
i,j = 1
The sequence in Table III represents different perturbation amplitudes (e) and scale contents (MI, Mu). We can identify a modon diameter with a perturbation wavelength L / M for wavenumber M = 5. Thus, there are perturbations with scales as large as or larger than the modon scale (expts. 10--14) or with scales smaller (expt. 17) or with scales spanning the modon scale (expts. 15, 16). The perturbation experiments have a weak dissipation (equivalent to that of expt. 9). The value of K2 is small enough to alter only slightly the modon propagation and persistence (in expt. 9 max ~"and Xc are within 10% of their inviscid values after At ~ 5), yet it provides a sink for enstrophy which can accumulate on the grid scale when e is large enough that a twodimensional turbulent cascade occurs. For weak perturbations (small e) modons are unaffected in their behavior, at least for the integration times tested. For example, the amplitude and central position of the modon in expt. 10 are the same at t = 5 as they are in expt. 9, to within an uncertainty no larger than the perturbation amplitude.
i
i
\ s
\
\
D
D
\ \
~struction~ \ \
¢
S
I
\
\ \ \
D \
~urvival~
Ol
o
\ S
r i
i 2
P 3
5 (M~MI) -I/2
Fig. 10. A regime diagram for the survival (S) or d e s t r u c t i o n ( D ) o f a barotropic m o d o n in the presence o f a r a n d o m perturbation. This diagram is based u p o n the results o f expts. 1 0 - - 1 7 . e is the area average RMS v o r t i c i t y a m p l i t u d e o f the perturbation ( t o be c o m pared w i t h a m o d o n v o r t i c i t y m a x i m u m o f 1 6 . 1 ) , and the perturbation scale is a geometric m e a n scale ( n o r m a l i z e d so that u n i t y is the m o d o n scale).
234 M=, M= - 1,5
rms = 0.37
MI,M= - 6 , 2 5
C . I . • 0.1 M,, M, • I, 25
rm$ , 0 . 0 1
rms • 0.08
C.1. • 0 . 0 3
C.1. • 0.01
o_
Fig. 11. Perturbation streamfunetion ~ ( x , y) for three sets of wavenumber intervals (Mu, M1) obtained by solving V2~ = Z, where Z is determined from eq. 8. For each panel the (area average) RMS value for • is listed.
For sufficiently strong perturbations, however, the m o d o n is destroyed as a separate entity by participation in a turbulent cascade. We have used a subjective criterion for m o d o n survival or destruction, y e t in our solutions the issue does n o t appear to be a subtle one because the modon, while surviving, is readily identified (the first panels in Figs. 12, 13). The critical perturbation amplitude required for destruction is a function
235 tmO
1~ 0.4
t =0.8
y
3(
Fig. 12. Vortieity patterns for an example of rapid modon destruction (expt. 14). The CI is 2.0.
of the scale c o n t e n t of the perturbation. Figure 10 is a regime diagram in p e r t u r b a t i o n vorticity amplitude and geometric mean horizontal scale. The critical amplitude decreases as the scale increases. Stated alternatively, vorticity on larger scales is more efficient in destroying m o d o n s than vorticity on smaller scales. Perturbation s t r e a m f u n c t i o n patterns f or e = 1 and various scales are shown in Fig. 11. The RMS amplitude o f • increases with pert urbat i on scale
236
tzl
t :10
Fig. 13. Vorticity patterns for an example of slow modon destruction (expt. 12}. The CI is 1.0.
237 approximately as its square, as expected from a second integral of vorticity. Thus, the boundary between m o d o n survival and destruction is less scaledependent when the perturbation amplitude is measured in terms of stream function rather than vorticity. This boundary occurs at an RMS ( e ~ ) value around 10--20% of the m o d o n e x t r e m u m in streamfunction. The mechanism for m o d o n destruction by perturbation appears to be advective shearing o u t of the m o d o n vorticity contours (which are material lines approximately) b y the large scale perturbations. Identifying the latter with the "energetic eddies", this mechanism is known to be active in twodimensional turbulence (Kraichnan, 1971). This process is indicated b y Fig. 12, which is an example of rapid m o d o n destruction. For smaller values of e, slightly above the critical value, the destruction rate is lower, and the shearing o u t of m o d o n vorticity contours seems to be selectively occurring at the edges of the vorticity centers and subsequently eroding inward. An early consequence of this is a separation of the two m o d o n centers, so that they cease to interact strongly and can begin to move westward and meridionally and otherwise act as isolated vortices (McWilliams and Flierl, 1979). This is illustrated in Fig. 13, which is an example of slow m o d o n destruction. The behavior shown in Fig. 13, which is based on expt. 12, is qualitatively the same in expts. 13 and 16 as well. Thus, the nature of the m o d o n destruction process in the slightly supercritical regime seems moderately independent of either scale content or the particular rand o m realization of the perturbation field. 7. CONCLUSIONS Numerical solutions of the barotropic, ~-plane potential vorticity equation have been examined to address the three aspects of m o d o n behavior listed in the introduction. In summary the results are as follows. (i) Modons can be successfully calculated b y standard numerical techniques if the resolution scales in space and time are sufficiently small. In particular, a b o u t 20 grid points per m o d o n diameter are required to obtain greater than 95% accuracy in the bulk propagation rate using second-order finite~lifference techniques. (ii) Under the influence of m o m e n t u m dissipation, modons decrease in amplitude, reduce their zonal propagation rate and expand their meridional scale. The first two processes occur in ways which are insensitive to gross aspects of the nature of the dissipation, and the third is simply related to the order af the dissipation law. After sufficient amplitude decay the m o d o n structures make a transition to a dispersive Rossby wave regime. While still within the m o d o n regime the decline in amplitude and speed crudely follows a m o d o n dispersion curve. (iii) Modon$ are resistant to perturbations of small amplitude and are destroyed b y perturbations of moderate amplitude. The critical amplitude for destruction is dependent u p o n the scale c o n t e n t of the perturbation, in a
238
manner consistent with advective shearing on larger than m o d o n scale being the dominant destructive mechanism.
ACKNOWLEDGMENTS
The numerical calculations were performed at the National Center for Atmospheric Research, which is sponsored by the National Science Foundation. Programming was done b y Dr. Julianna H.S. Chow. The international collaboration which this paper represents arose under the auspices of the POLYMODE Program.
REFERENCES Arakawa, A., 1966. Computational design for long term numerical integration of the equations of fluid motion: Two dimensional incompressible flow. Part I. J. Comput. Phys., 1 : 119--143. Berestov, A.L., 1979. Solitary Rossby waves, Izv. Acad. Sci., U.S.S.R., Atmos. Oceanic Phys., 15 : 648--654. Flierl, G.R., Larichev, V.D., McWilliams J.C. and Reznik, G.M., 1980. The dynamics of baroclinic and barotropic solitary eddies. Dyn. Atmos. Oceans, 5: 1--41. Kraichnan, R.H., 1971: Inertial-range transfer in two- and three-dimensional turbulence. J. Fluid Mech., 47: 525--535. Larichev, V. and Reznik, G., 1976. Two-dimensional Rossby soliton: an exact solution. Rep. U.S.S.R., Acad. Sci., 231 (5) (also POLYMODE News, 19). Lilly, D.K., 1965. On the computational stability of numerical solutions of time-dependent non-linear geophysical fluid dynamics problems. Mon. Weather Rev., 9 3 : 1 1 - - 2 6 . McWilliams, J.C., 1980. An application of equivalent modons to atmospheric blocking. Dyn. Atmos. Oceans, 5 : 43--66. McWilliams, J.C. and Flierl, G.R., 1979. On the evolution of isolated, non-linear vortices. J. Phys. Oceanogr., 9: 1155--1182. Orzsag, S.A. and Jayne L.W., 1974. Local errors of approximate solutions to hyperbolic problems. J. Comput. Phys., 14: 93--103. Scott, A.C., Chu, F.Y.F. and McLaughlin, D.W., 1973. The soliton : a new concept in applied science. Proc. IEEE, 61: 1443--1483. Stern, M.E., 1975. Minimal properties of planetary eddies. J. Mar. Res., 33: 1--13. Thompson, P.D., 1961. Numerical Weather Analysis and Prediction. Macmillan, New York, NY, 170 pp. Zabusky, N.J. and Kruskal, M.D., 1965. Interaction of "solitons" in a collisonless plasma and the occurrence of initial states. Phys. Rev. Lett., 15: 240--243.