Spatio-temporal stability analysis applied to monsoon anticyclone flow

Flow on a beta-plane driven by a steady localised anticyclonic forcing of potential vorticity (or equivalently a mass source) is considered as a simple model of the Asian monsoon flow in the upper troposphere. Previous authors have noted that the response may be steady, or unsteady, according to the magnitude of the forcing, with the unsteadiness manifested as westward eddy shedding. A detailed study of the transition between steady and eddy-shedding regimes reveals a third regime ('break up'), for intermediate forcing magnitude, where the flow is steady in the neighbourhood of the forcing, but the westward extending plume of low potential vorticity breaks up into isolated anticyclonic vortices some distance away from the forcing region. A related spatio-temporal instability problem for flow on a beta plane is specified and analysed. The flow can be stable, convectively unstable or absolutely unstable. It is argued that these three stability regimes correspond to the steady, break-up and eddy-shedding regimes for the forced flow and good quantitative correspondence between the regimes is demonstrated by explicit solution of the spatio-temporal stability problem.

The upper-level anticyclone has broader importance for several reasons, one being its remote influence on weather and climate, another being its role in transport from troposphere to stratosphere. In the lower part there is fast convective transport from the surface to the upper troposphere. Then from the upper troposphere air may be taken up into the tropical stratosphere via the large-scale upwelling circulation, or alternatively be transported laterally into the extratropical lower stratosphere through disturbances to the subtropical jet that depend on the time variation of the anticyclone itself. Thus the lower and upper parts of the monsoon combine to form a potentially rapid pathway for chemical species and aerosols originating from anthropogenic emissions at the surface in the south and east Asian regions to reach the stratosphere, with implications for climate and for ozone, since the chemical species include short-lived halogen species (Park et al., 2009;Randel et al., 2010;Bourassa et al., 2012;Fiehn et al., 2017).
Observational studies have revealed that an important aspect of the time-dependent behaviour of the upper-level anticyclone is what is described as "vortex shedding" or "eddy shedding". The upper-level forcing provided by the latent heating generates a region of low PV. Rather than reaching some equilibrium size, or simply increasing in size, this region is seen to break up, forming an anticyclonic eddy which moves away from the forcing region most commonly to the west ("westward shedding"), but sometimes to the east (e.g., Vogel et al., 2015). The focus in this paper will be on the westward shedding.
Using re-analysis data, Hsu and Plumb (2000) illustrated a particular westward-shedding event occurring during the monsoon season 1990, visible as coherent PV structure breaking off the bulk PV-low associated with the main monsoon anticyclone. Popovic and Plumb (2001) showed that shedding events are not rare, but typically happen several times during a monsoon cycle. They further found the shed eddies to be rather shallow and confined to the upper troposphere/lower stratosphere (UTLS). More recently Garny and Randel (2013) examined the variability of the anticyclone, including the PV structure, and its association with variability in carbon monoxide.
The observed variations (typically on a time-scale of about two weeks) in the monsoon anticyclone have been interpreted in different ways. Zhang et al. (2002) have suggested a "bi-modality" with eastward or westward displacements of the centre of the anticyclone. Vogel et al. (2015) have described a "split-state" where there are two pronounced and zonally separated PV extrema. Nützel et al. (2016) investigated multiple re-analysis datasets and found that only one (NCEP-1) out of seven showed pronounced and consistent signs of bi-modality. It seems quite possible that the underlying process is eddy shedding, which may be interpreted as bimodality or split-state according to the particular diagnostic and perhaps dataset considered.
The potential importance of the spatial and temporal variability of the monsoon anticyclone for chemical transport and the meteorological interest in the phenomenology of the variability implies a need for dynamical understanding. One basic question, considered by Garny and Randel (2013), for example, is the extent to which the temporal variability of the monsoon anticyclone is directly determined by the variability of precipitation and the associated latent heating versus the extent to which it is an intrinsic feature of the dynamics of the upper-tropospheric flow.
Some aspects of the structure of the monsoon circulation, both at lower and upper levels, are captured by the linear dynamical model of Gill (1980), but because this model relies on frictional damping, which is difficult to justify in the upper troposphere, its relevance to the structure of the upper-level anticyclone is questionable (e.g., Sardeshmukh and Hoskins, 1985;Lin et al., 2008). Furthermore the model is steady and therefore its relevance to any time dependence of the anticyclone is likely limited. A different approach has been taken by Hsu and Plumb (2000), who consider a single layer shallow-water model, representing the upper troposphere, with a specified localised anticyclonic PV forcing (or equivalent a localised mass source). Hsu and Plumb (2000) showed in numerical simulations that the phenomenon of westward eddy shedding can spontaneously emerge as response to a steady localised forcing and hence provide a possible mechanism for time dependence that does not rely on the variability of the forcing. They further investigated the transition between a steady state of the system and a time-periodic eddy-shedding state when imposing a meridional background PV gradient for a given forcing. This transition, however, is fundamentally different from the transition investigated in this paper, since Hsu and Plumb (2000) find steady behaviour in the case of strong forcing magnitudes (or equivalently weak background PV gradient) and unsteady behaviour for weak forcing magnitudes.
A very similar fluid dynamical model, but motivated by oceanographic problems, for example the formation of isolated eddies from the Mediterranean outflow, had previously been studied by Davey and Killworth (1989, hereafter DK89). They found the response to transition from a steady, zonally symmetric response extending westward from the forcing region (a so-called beta-plume) at small forcing amplitude to a westward-extending series of periodically shed vortices at larger forcing amplitude. A similar behaviour was found by Cenedese and Linden (1999) in a series of laboratory experiments to study the response to a continuous mass flux into a two-layer system with a sloped bottom topography (inducing a beta-effect). They also found a transition from a zonally symmetric response to periodic shedding when increasing the flux. The phenomenon of eddy shedding from a region where there is a localised PV forcing is therefore potentially important to a wide range of problems in geophysical fluid dynamics.
The numerical and laboratory experiments in the studies by DK89, Cenedese and Linden (1999) and Hsu and Plumb (2000) were all performed on a domain with relatively small zonal scale. This restriction on the domain size may seem consistent with the spatial scales of the problem (motivated by the monsoon anticyclone flow), but can make it difficult to observe the full range of behaviours within the system since in certain situations important features of the forced response only develop a large distance to the west of the forcing region (Section 3).
The study reported in this paper will analyse the mechanisms controlling westward eddy shedding from a steady, localised anticyclonic forcing in terms of a spatio-temporal stability analysis. The concept of spatio-temporal instabilities of shear flows was first introduced in the field of plasma physics (e.g., by Briggs, 1964;Bers, 1973) and finds application in many branches of fluid dynamics, particularly in the evolution of spatially evolving shear flows and their response to excitation (Huerre and Monkewitz, 1990). The possible application to a range of problems in geophysical fluid dynamics was quickly recognised (e.g., by Merkine, 1977). Pierrehumbert (1986) applied the idea of spatio-temporal stability to the Charney problem and investigated the formation of baroclinic instability for different surface wind magnitudes. Lin and Pierrehumbert (1993) extended this approach to flows with both latitudinal and vertical structure. Schär and Smith (1993) performed a spatio-temporal stability analysis to link eddy-shedding phenomena caused by a flow over an isolated mountain in a non-rotating shallow-water model to the absolute instability of certain idealised steady solutions. In a meteorological setting, Hohenegger et al. (2006) argued that the lack of predictability skill for certain precipitation events in numerical models is not fully explained by problems in modelling moist convection, but is partially caused by a spatio-temporal amplification of small-amplitude perturbations. Wang (2011) performed a comprehensive spatio-temporal stability analysis to explain small-scale structures observed in oceanic eastern boundary currents. Diaz and Aiyyer (2015) used the idea of spatio-temporal instabilities to model large-scale wave patterns found within the African easterly jet.
The structure of the paper is as follows. The specifics of the model used in the present study is described in Section 2. In Section 3 we will discuss a range of behaviours observed in steadily forced numerical experiments. This will in particular involve a transition (depending on details of the forcing) between two distinct states of the system in which the response to the steady forcing spontaneously develops a temporal dependence. As we will show, the transition is related to a change in stability properties of the system. We will therefore subsequently analyse the stability of an idealised representation of the forced flow in two different ways. In Section 4 we perform a spatio-temporal stability analysis based on the Rayleigh equation of the system. This is a delicate calculation and therefore in Section 5 we present a complementary analysis of the stability of the flow using direct numerical simulations and compare our findings to the theoretical predictions obtained in the previous section. Section 6 will then use these theoretical predictions to explain aspects of the behaviour of the response to a localised steady forcing. Finally we will summarise and discuss our main findings in Section 7.

MODEL
In the following we will consider a single-layer beta-channel model, intended as a simple representation of the subtropical Asian monsoon system. A single-layer model of the low-latitude upper troposphere has been used in several previous studies (e.g., Hsu and Plumb, 2000;Kraucunas and Hartmann, 2007;Bao and Hartmann, 2014) and can be justified on the basis of the "tropical scaling" first set out by Charney (e.g., Vallis, 2017, chapter 18.8) which predicts that, in the absence of diabatic heating, the motion of each horizontal layer of fluid is essentially independent of those above and below. A restriction relative to previous studies is that, rather than considering the full single-layer equations, we make the quasi-geostrophic (QG) approximation. This is justified on the basis that the main effect of the approximation is to omit gravity waves, including equatorial Kelvin waves, and that the subtropical vortex-shedding phenomenon does not depend significantly on such waves. However, note that Hsu and Plumb (2000) consider some tropical cases where there is a more substantial role for Kelvin waves. Relaxing the QG approximation will of course alter some of the quantitative details of the behaviour, but it should also be noted that the single-layer model representation is itself approximate and will not correspond precisely to the behaviour of a three-dimensional atmosphere. Diabatic heating is included in the model only to the extent that it provides a specified anticyclonic forcing since we are interested in the flow response rather than the details of the forcing itself.
The non-dimensional model equations are given by: where q describes (QG) PV, is the stream function, u = − y and v = x are the fluid velocities in the zonal (x) and meridional (y) directions, respectively, F is a PV forcing (or equivalently a mass source), S is a space-dependent damping, is a hyperdiffusion parameter and ∇ is the two-dimensional gradient operator. The parameter̂represents the non-dimensional meridional gradient of background PV. Equation (1) is non-dimensionalised, scaling horizontal lengths by the Rossby radius of deformation R D , times on a scale of 1 day, and all other quantities correspondingly.
For the purpose of the results to be presented, we will takê= 2.2, which corresponds to about 20 • latitude, when assuming R D = 1, 250 km. This value of R D is further consistent with dynamical structures of the Rossby height of roughly H R = fR D N −1 ≈ 3 km, assuming f ≈ 0.5 ⋅ 10 −4 s −1 and N ≈ 2 ⋅ 10 −2 s −1 , which are values appropriate for the subtropical upper troposphere. In the set of single-layer equations, this further leads to a layer thickness of about H = R 2 D f 2 ∕g ≈ 400 m, with g being the gravitational acceleration. (For comparison, Kraucunas andHartmann, 2007 used a value of H = 200 m andBao andHartmann, 2014 used H = 600 m.) The model domain is taken to be −240 < x < 80 and −20 < y < 20. Given the non-dimensionalisation described above, this domain is clearly much larger than can be justified by application to the subtropical atmosphere. However, the objective here is to remove, as far as possible, the effect of rigid boundaries in y and the effect of the periodic domain in x, to allow a clear understanding of the behaviour of the system. Note that most figures will show the response only within a part of the full domain. Equation (1) is numerically integrated using a spectral representation in x, truncated at wave number 2,048, and a second-order central-difference scheme with 256 equally spaced grid points in y. We further assume = 0 at the meridional boundaries, implying vanishing normal flow. In time the model is integrated via a second-order Runge-Kutta method with 200 steps per day. To ensure numerical stability of the model, was chosen to damp the shortest waves resolved on a time-scale of about 100 days. The spectral representation in x implies periodicity, but we wish to consider the solutions to Equation (1) on an infinite x domain, with no communication between large positive x and large negative x. To achieve this, S(x) is set to be non-zero only in a sponge region with x < −200, which is to the far west of the forcing region (or the east, due to the periodic boundary conditions). Note that, for the steadily forced experiments in Section 3, the damping acts on the full PV field, while for the experiments in Section 5 it acts only on the zonally varying part of the response. The sensitivity of the solutions to Equation (1) to the precise details (structure and magnitude) of the sponge region S(x), the hyperdiffusion parameter , the domain size and the spatial resolution was thoroughly tested and all the results presented below are robust to changes in these various quantities.
The forcing is given by a confined mass source with spatio-temporal structure F(x, y, t), switched on smoothly over the first 10 days of each experiment and defined via Here, F 0 controls the forcing magnitude, r 0 the non-dimensional length-scale of the forcing and 0 the time dependence. The primary motivation for the work reported in this paper is to understand the behaviour of the system starting from an initial condition at rest, when the forcing is steady ( 0 = 0). But in interpreting this behaviour, it is useful also to consider the response to a time-dependent forcing ( 0 ≠ 0), starting from an initial condition of a specified zonally uniform flow. In the case with 0 = 0, the forcing represents a constant mass source into the system, but note that the dynamics described by Equation (1) stay invariant when adding a constant to and thus a statistical steady state can be reached without the explicit need for dissipative processes. Further details on the nature of the forcing will be given in Sections 3 and 5.

DIFFERENT TYPES OF RESPONSE TO A STEADY FORCING
As noted in Section 1, the horizontal structure and temporal evolution of the response to a continuous ( 0 = 0 in Equation (2)) mass source in a single-layer model has previously been investigated by DK89, who observed a state transition when increasing the forcing magnitude while keeping the length-scale of the forcing fixed. They showed that the response changes from a steady, westward extended plume of PV anomaly to a state in which vortices form spontaneously and periodically inside the forcing F I G U R E 1 Day 400 PV snapshots (shading) using different combinations of forcing magnitudes F 0 and length-scales r 0 , labelled cases (a)-(f). The system is exhibiting various states of temporal behaviour and spatial structure. We can distinguish steady, break-up and shedding states. The contours show the positive part of the stream function field, using different contour intervals for the different plots. The forcing used in (b), (c), (e), and (f) satisfies the shedding criterion in Equation (A3), while the forcing in (a) and (d) does not. The superimposed ellipses indicate the forcing region region and are subsequently shed to the west, forming a zonal chain of discrete eddies.
Furthermore DK89 gave a theoretical prediction of the threshold value of the forcing at which the state transition takes place, based on consideration of whether the PV gradient resulting from the steady linear response became negative, hence satisfying the Rayleigh inflexion point criterion necessary for instability. They argued that, once this criterion was satisfied, the steady response would no longer be observed and instability would result in the formation of discrete eddies. A derivation of an equivalent criterion for our QG model and a specified forcing profile is given in the Appendix.
We now investigate whether the behaviour identified by DK89 and the shedding criterion (Equation (A3) in the Appendix) applies in our system and how it depends on the parameters of the forcing. We consider Equation (1) integrated forward in time from an initial state at rest. Figure 1 shows the model response at day 400 for the specified forcing profile and different sets of forcing parameters, which we will refer to as cases (a) to (f) within this section. The response for two different forcing length-scales r 0 are shown, using three different values of forcing magnitude F 0 each (Equation (2)). Note that DK89 used a fixed forcing radius of r 0 = 3, that is, a forcing radius that is somewhat larger than the deformation radius, as in the examples presented here, but that their value of̂and various details of the underlying model were different. However, our numerical simulations show a good match with DK89 when the corresponding physical parameters are chosen appropriately. This provides reassuring evidence that the QG approximation, which was not used by DK89, does not disrupt the eddy-shedding behaviour.
For weak forcing magnitudes (cases (a) and (d)) the response is steady and does not show any sign of zonal variability west of the forcing. The corresponding forcing profiles do not satisfy Criterion A3 and the response does not have a reversed mean meridional PV gradient, as can be seen in Figure 2. The flow in these cases is, as DK89 have shown, well described by steady, linear theory.
For forcing profiles that satisfy Criterion A3 (cases (b), (c), (e) and (f)), we observe the formation of discrete vortices despite the use of a steady mass source, as described and predicted by DK89. Figure 2 shows that the corresponding cases are associated with a region of negative (or vanishingly small in case (b)) mean meridional PV gradient. Note that cases (e) and (f) look very similar in terms of the mean PV gradient, but the detailed structure seen in Figure 1 is very different.
It is here important to remember that the Davey-Killworth criterion in terms of a sign change in PV gradient is simply a necessary condition for instability that is entirely based on a steady solution to the problem. The response in cases (b), (c), (d) and (f), on the other hand, is unsteady and thus the significance of a sign change in the mean PV gradient is not so clear. The profiles displayed in Figure 2 show a time average which accounts for times when there is an eddy at x = −10, and times when there is not. However, the idea that the time-mean PV gradient changes signs in the non-steady cases is consistent with the idea of the forcing generating a PV plume that is subsequently breaking up into discrete vortices due to the growth of small random perturbations because of some instability.
In contrast to the findings of DK89, we observe two qualitatively different kinds of behaviour when the forcing profile satisfies the shedding criterion given in the Appendix (Equation (A3)). For certain parameter combinations (case (e) in Figure 1) we find the response in and close to the forcing region to be almost steady and qualitatively similar to the solution predicted by steady, linear theory, but with a sustained meridional PV gradient reversal, as shown in Figure 2. In these cases the westward extending plume becomes perturbed and breaks up F I G U R E 2 Time mean meridional PV gradients along x = −10 for the same experiments as in Figure 1, for length-scale r 0 = (a) 2 and (b) 4, and with various values of forcing magnitude F 0 (a) (b) into discrete patches a finite distance away from the forcing. The corresponding state is therefore characterised by weak time dependence near the forcing and strong time dependence to the far west of the forcing. Note that such behaviour might not be observed for small domain sizes or short time integrations where the flow does not have the opportunity to extend sufficiently far westward and might incorrectly be classified as completely steady. In other cases (e.g., case (f)), the vortices form directly within the forcing region and are then shed westward. This state is characterised by a strong time dependence of the flow near the forcing region (as well as further to the west). In the following, we will refer to the three different types of system behaviour as steady, break-up and shedding states.
To identify the type of response to a certain forcing with given scale r 0 and magnitude F 0 , we evaluated the temporal behaviour close to the forcing (at x = −10) and the zonal structure of its stream function at day 300. In the steady cases we expect no strong variations, neither in time nor in space. We classify the system as being in the break-up state if the stream function stays steady near the forcing region but develops strong zonal and temporal variations far to the west. The spatial amplification of perturbations in this state suggests the system to be convectively unstable, as explained in detail in Section 4. Finally we classify the response as being in the eddy-shedding regime if we observe strong temporal perturbations near the forcing and a consistently perturbed zonal structure west of the forcing. We will later in Section 6 use the characteristics of the spatial structure and temporal variability that have been discussed in this paragraph to classify the state of the response.
The evolution of the system in the shedding state typically leads to the formation of an almost periodic chain of vortices, as can be seen in Figure 1b,c,f. Note a significant meridional drift of the formed discrete vortices on top of the predominantly zonal propagation, which is particularly strong for small and strong eddies. This drift is a result of the beta-effect, as discussed, for example, by Davey and Killworth (1984). To reduce the impact of this meridional drift on our diagnostics, many of the results shown will be meridionally integrated over a range of roughly the meridional extent of the forcing region. The shed vortices seem to travel with a relatively constant zonal velocity, as also described by other authors (Davey and Killworth, 1989;Hsu and Plumb, 2000). The periodic formation of vortices therefore leads to a periodic signal of the response at fixed locations west of the forcing. To further quantitatively characterise the shedding process, we estimated the frequency at which vortices are created and shed from the forcing region (the shedding frequency) by computing the temporal Fourier power spectra of the stream function response to the west of the forcing for experiments (b), (c), (e) and (f) shown in Figure 1. Although the different experiments cover a wide range of parameter combinations for the forcing, and in addition include responses in the shedding and break-up states, the various power spectra (not shown) all indicate dominant frequencies in a relatively narrow range of about 0.8 < < 1.1.
The structure of the response in the break-up state shown in Figure 1e suggests a zonally amplifying perturbation on top of an almost steady beta-plume. We therefore propose that certain aspects of the response to a localised steady forcing can be explained as the the result of a spatio-temporal instability of the forced beta plume that is itself well described by steady linear theory (DK89; see also the Appendix to this paper). In particular, we will show that the transition between the steady, break-up and shedding states can potentially be explained by a change in stability properties of the system. The interpretation of the shedding process as a result of a spatio-temporal instability further potentially allows us to predict the dominant frequencies that we reported in the previous paragraph.

STABILITY OF CONFINED PARALLEL SHEAR FLOWS
We will in this section investigate to what extent the three qualitatively different states identified in Section 3 can be explained as the result of changes in the spatio-temporal stability properties of the system. We will therefore perform a comprehensive spatio-temporal stability analysis for a set of basic flows representing the steady parallel shear flow that is predicted as the beta-plane steady linear response to a localised vorticity forcing far to the west of the region where the forcing is applied. The QG PV model described by Equation (1) can be linearised for a zonal wave-like perturbation around the zonally uniform basic state withū(y) = − y (y), that is, where ′ describes the meridional structure of the wave, k is the x-wave number and the frequency. In principle, both k and can be complex. The corresponding complex phase speed is then given by c = ∕k. The standard approach to analysing stability is to write the resulting linearised equation in the Rayleigh form together with boundary conditions ′ = 0 in y. Equation (3) can be used to define a function D(k, ), where the condition D = 0 represents the dispersion relation of the system and links the (complex) wave number and frequency of a supported wave. Note that D(k, ) = 0 generally has multiple solutions, each corresponding to a different wave mode supported by the system. Determining D(k. ), and thus determining the spatio-temporal stability properties, will require numerical solution of Equation (3), which we describe in Subsections 4.2-4.4 below.

Theoretical background
In this subsection we summarise the relevant theoretical background on temporal, convective and absolute instability. The stability behaviour is determined by D(k, ) = 0.
Since both k and can in principle be complex, we can distinguish between different types of instability. The familiar case of temporal instability corresponds to real k and there is instability if the corresponding has a positive imaginary part. Another possibility is what is known as convective instability, where a mode grows in space but not in time, that is, k has a non-vanishing imaginary part and is real. Note that positive k implies growth in the negative x direction and negative k growth in the positive x direction.
It is now understood that the response of a flow to a localised forcing is not completely determined by whether or not there is temporal or convective instability. Briggs (1964) and Huerre and Monkewitz (1990), for example, determined this by considering the evolution of a localised wave packet. If there is temporal instability, the wave amplitude may grow, but if the wave packet propagates away from its initial location sufficiently fast, then the waves will be observed to decay in amplitude at the initial location. This situation corresponds to convective instability. The alternative is that the waves at the initial location are observed to increase in amplitude for large times. This situation is called absolute instability and corresponds to both k and having non-vanishing imaginary parts. The wave packet interpretation makes it clear that absolute instability has to be associated with a saddle point in the complex dispersion relation at an isolated frequency a and corresponding wave number k a , since such a point corresponds to a wave mode with vanishing group velocity d ∕dk = 0 which allows for energy to grow in place.
The distinction between convective and absolute instability is essential to predict the behaviour of disturbances to a flow that will be observed in some finite range of x. Determining where there is absolute instability requires a detailed understanding of the behaviour of the roots of D(k, ) = 0 in the complex plane. Note that the interpretation of convective and absolute instability in terms of an initial wave packet perturbation still holds when considering a problem where a basic flow is perturbed by a spatially localised periodic forcing with fixed frequency, which corresponds to the idea of a continuous formation of wave packets with modulated amplitude at a fixed location. Such an interpretation is conceptually closer to the problems considered in this paper and forms the motivation for the experiments in Section 5.
It is not sufficient to analyse the form of D locally in the (k, )-plane to determine whether a flow is convectively or absolutely unstable. As discussed by Briggs (1964), the convective roots of the dispersion relation, that is, roots with real and k complex, only correspond to a physically relevant disturbance if k( ) crosses the real k-axis as decreases from large positive values towards zero. It follows as a corollary that there is convective or absolute instability only if there is temporal instability. As further discussed by Briggs (1964), there is a similar criterion for absolute instability to determine whether or not the corresponding saddle point (which is formed as a double root of D) corresponds to a physically relevant disturbance. The condition is, again, based on the global behaviour of the dispersion relation D around the saddle point at k a = k( a ) and is usually described in terms of curves in the complex k-plane defined by the (inverse) dispersion relation k( ) with i = fixed and r = ℜ varying, often called Bromwich images. The corresponding lines in the complex plane with fixed i are called Bromwich contours. There may be several Bromwich images for a given Bromwich contour (i.e., a specified i ) if k( ) is multivalued. A saddle point (or double root) of D now only describes a physically relevant mode of absolute instability if it is formed by two merging Bromwich images that originate from opposite half-planes in k-space as i is varied, as was particularly stressed by Pierrehumbert (1986). It is clear from the above that finding convectively (and absolutely) unstable modes for the Rayleigh equation of our system (Equation (3)) and similar equations is a subtle undertaking and different approaches to this have been developed and used. Some authors (e.g., Pierrehumbert, 1986;Polvani and Pedlosky, 1988;DelSole, 1997) were able to formulate a problem for which the dispersion relation could be obtained in analytic form. Lin and Pierrehumbert (1993) considered the stability of a zonally symmetric jet in a QG beta-plane model which leads to an eigenvalue problem in k and c, similar to Equation (3), and pointed out the subtlety of these kinds of problems if the dispersion relation cannot be obtained in analytic form. They described an iterative procedure starting by obtaining the temporal dispersion relation (found by imposing real k in Equation (3)) and subsequent tracking of the unstable mode as k becomes complex. This method is conceptually similar to the approach we will be using, described below, however, it can only be used to identify absolutely unstable modes and does not directly give any information on pure convective instability of the system or other characteristics of the dispersion relation (in contrast to our approach).

Calculation of temporal stability
In the rest of this section we will progress towards a full solution of the spatio-temporal stability problem to determine whether, for a specified basic state, there is absolute instability. This is practically done by solving Equation (3) which, together with the boundary conditions, can be regarded as a generalised eigenvalue problem that straightforwardly implies c given k (and hence given k), or as a different eigenvalue problem that straightforwardly implies k given c = ∕k. However it is not so straightforward to determine k given . We solve Equation (3) numerically by simple discretisation using second-order centred differences and then using standard MATLAB routines to solve the resulting matrix eigenvalue problem. The eigenvalue of interest has to be selected from the set of all eigenvalues of the matrix. The criteria used for this selection are described later in this subsection. All the results reported use a discretisation with 256 points equally spaced over a domain of size 40 in the y-direction.
Robustness of the results to changes in the number of points and the size of the domain was checked.
For illustration of our methods, we consider the particular family of basic state profiles defined by the stream function for |y| < , where determines the strength of the flow and its meridional width. The basic state is supposed to model a zonally extended steady plume, motivated by the steady linear solution of Equation (1) (as described in Section 3). The form of in Equation (4) has purposely not been chosen to match the steady linear response to the forcing specified by Equation (2), since that would imply a discontinuity in the corresponding PV field at the edges of the forcing region and hence numerical difficulties. However, our aim is to explain the response to a steady localised forcing, as described in Section 3. In practice in the numerical simulations nonlinear effects seem to smooth any discontinuities and the smoothed discontinuities then manifest as secondary peaks in the meridional PV gradient at the edges of the forcing region ( Figure 2). Hence profile 4 turns out to be a useful semi-quantitative representation of the profiles seen in the simulations, as can also be seen from Figure 3 and the corresponding explaining paragraphs in Section 6.
The general structure of the basic state is illustrated in Figure 3, also showing the background zonal windū(y) = − y and meridional PV gradient q y =̂+ū − 2 yū . Note that the meridional PV gradient is single-signed for < G, where G( ,̂) can in principle be calculated exactly from the basic state profile in Equation (4) and the definition of PV in Equation (1). In cases where < G, we expect the system to be stable in a temporal sense, which further implies convective and absolute stability. For a fixed =6 this is the case for roughly <5.
Since the analysis is relatively complicated, we begin by considering the temporal stability problem, which is solved by specifying a (real) value for k and computing, via the generalised eigenvalue problem in Equation (3), the (complex) phase speed c, which immediately implies the corresponding (complex) frequency = ck. The fastest growing mode is that, considering all (real) values of k, with the largest positive value of i = . For the basic state shown in Figure 3 with specified choices of the parameters and , it was straightforward to solve the eigenvalue problem for c(k) numerically using the method outlined above for many choices for imposed, real k. As noted previously, a solution of the matrix eigenvalue problem for a given k gives many eigenvalues (the same F I G U R E 3 (a) Idealised basic state profile with = 9.5 and = 5.4 defined via Equation (4) and associated zonal wind u = − y and meridional PV gradient q y =̂+ū − 2 yū . (b) measured state profile extracted from a direct numerical simulation with forcing parameters r 0 = 4, F 0 = 1.3 (defined later in Section 6) number as the number of points used in the discretisation). The selection procedure applied, for each k, was to choose the eigenvalue with largest imaginary part (most of the eigenvalues had zero imaginary part and represent an approximation to the "continuous" spectrum of the Rayleigh problem). This specifies a unique c(k). Figure 4 shows the temporal dispersion relation in terms of frequency and wave speed c as a function of (real) k for different combinations of basic state parameters and . We find the system to be temporally unstable for a range of wave numbers in all three cases shown. The growth rate i seems to generally increase with increasing , specifying the strength of the basic state flow, and decreasing , specifying the width of the flow. The real part of the frequency shows relatively little dependence on either or and an almost inverse linear dependence on k gives rise to a relatively constant value for c r = ℜc in all three cases.
A distinct, continuous dispersion relation can be identified in all cases where the PV gradient changes sign, which we interpret as a coherent and physically relevant temporal mode. The overall shape of the structure is relatively similar for the various parameter combinations and shows a clear local maximum in i . However, the value of the maximum growth rate varies and seems to increase for narrower and stronger basic states, making the flow less stable in these cases. Further note that the range of k with i > 0 seems to decrease significantly as increases or decreases (consistent with the idea that the system becomes stable for large and/or small , as explained earlier). As noted previously, temporal instability of the system is a necessary condition for convective or absolute instability. In Section 4.3 we will use this condition to solve the convective instability problem starting from an unstable solution of the temporal problem. Figure 4 can be used to identify suitable temporally unstable modes, that is, a set of real wave number and corresponding complex frequency, that can be used to calculate potential convectively (or absolutely) unstable modes. Note that solving the eigenvalue problem in Equation (3) for imposed real k together with the selection process described above can lead to problems due to numerical uncertainties in cases where the corresponding solutions (i.e., the calculated phase speeds and thus frequencies) have small imaginary parts. The profiles shown in Figure 4 therefore have to be interpreted with caution in regions with small i and are most reliable around the corresponding local maximum of i (k). For the analysis to be presented in Sections 4.3 and 4.4, we therefore chose to start the calculations with values of k close to the local maximum in i .

Calculation of convective instability
Convective instability can imply growth in the positive or negative x-direction. For simplicity (bearing in mind the results shown in Figure 1) we shall here consider only F I G U R E 4 Temporal dispersion relation in frequency and wave speed for idealised basic states with different widths and magnitudes growth in the negative x-direction. This requires a solution to Equation (3) with positive imaginary wave number k i = k, satisfying the additional condition that it corresponds to a root that crosses the real k-axis as i decreases from large positive values towards zero. Both conditions require knowledge of the structure of the dispersion relation for varying . However, as noted previously, Equation (3), or its numerical representation, does not straightforwardly allow us to determine k( ), but rather to determine k(c) and then compute = kc, which increases the difficulty of finding physically relevant convectively unstable modes.
The two conditions mentioned above suggest the following approach to finding the required complex k( ) for real . We start with solutions to the temporal instability problem, each providing a real k(c) for particular c, with the corresponding i greater than zero. We then vary c in such a way to reduce i to zero. This can be done using the linearisation which allows, for example, Δc to be chosen such that Δ has zero real part. This gives us a constructive relation between the real and imaginary parts of Δc leading to the desired behaviour of Δ : where M = c(Δk∕Δc) + k and M r and M i are the respective real and imaginary parts. Note that, due to the assumption of linearity, we can change the sign of the real part of the change in frequency Δ i by changing the sign of the corresponding change in phase speed Δc. This allows us to specify, for example, a Δ with vanishing real part and negative imaginary part. Following similar considerations as for Equation (6), we can construct Δc in a way that leads to a change in frequency Δ with vanishing imaginary part: Starting from a known solution of the temporal problem (real k and with positive imaginary part), we can use Equation (6) to gradually and iteratively decrease the imaginary part of until it becomes real. While doing so, we track the corresponding value of k in the complex plane. This way we can construct a solution to Equation (3) with complex k and real that, if the corresponding k i > 0, corresponds to a convectively unstable mode.
Once a convectively unstable mode was identified, we can use Equation (7) to obtain the convective dispersion  (6) and (7) and a suitable solution of the temporal problem. The triangles show frequency/wave number pairs deduced from the numerical simulations discussed in Section 5 and provide a useful check on the validity of the approach to solving the dispersion relation relation by varying Δc in such way that the corresponding i stays zero, but r (and thus k) changes. Figure 5 shows the convective dispersion relation found following the described iterative method for three basic states with fixed width = 6 and different magnitudes .
We first consider the case = 6.5, for which Figure 5 shows a distinct local maximum in convective growth rate k i at r ≈ −0.85 with the corresponding real part of the wave number being k r ≈ 0.6. Following previous work on spatial instability in other fluid dynamical problems (e.g., Ho and Huerre, 1984), we interpret the convective dispersion relation as providing information on the growth of random perturbations to a zonally extended and steady beta-plume which results in the break-up behaviour identified in Figure 1e. We envisage that these perturbations correspond to a range of frequencies and that the perturbations that dominate as x decreases (i.e., as distance to the west of the forcing region increases) are those with largest spatial growth rate k i . Alongside the spatial growth the perturbations have an oscillatory structure in x with scale predicted by k r , corresponding to the x-scale of the eddies seen in Figure 1e, at least in the region where the evolution is well described by linear theory. Note that, as discussed in Ho and Huerre (1984) and similar papers, the convective dispersion relation also provides information on the forced problem, where r is imposed, for example, in a laboratory flow by adding an oscillatory forcing, but that is not relevant to the case we are considering here.
We investigated the sensitivity of the frequency and wave number of the convectively unstable mode with maximum spatial growth rate and found that, for varying basic state strength parameter 5 < < 7, their real parts stay relatively constant around values of r ≈ 0.8 and k r ≈ 0.6. The imaginary part of the wave number, however, seems to be zero for ≈ 5 and gradually increases with taking a value of k i ≈ 0.1 at ≈ 7. Recall that for a fixed = 6, as explained earlier, the flow is stable only when the meridional PV gradient is single-signed, which is the case for roughly < 5, which is in agreement with a vanishing convective growth rate as → 5. Now considering the cases = 8.6 and 13, it is seen from Figure 5 that k i ( r ) does not have an apparent local maximum, at least in the range of r considered. In fact it will be shown in Section 4.4 that for these values of there is absolute instability. As has been noted previously, absolute instability is associated with a saddle point in the dispersion relation D(k, ) with the transition from stability to instability as this saddle point crosses the real -axis into the positive i -half-plane. Since the saddle point is the result of a merging of two branches of the dispersion relation originating from opposite sides of the real k-axis (discussed in detail in Section 4.4) the iterative method for convective instability outlined above, following a certain branch of D with i = 0, can start to follow different branches of the function k( ) when the real frequency exceeds the real part of the frequency of the absolutely unstable mode (i.e., for | | > a,r ≈ 0.8, as explained in Section 4.4). When the basic flow is absolutely unstable, we therefore cannot be certain that the calculation of k i , that is, the results shown in Figure 5, is correct for large frequencies r . However, in these cases the flow behaviour at large times will be dominated by the absolutely unstable modes and the properties of the convectively unstable modes will not be of interest.

Calculation of absolute instability
To compute absolutely unstable modes, we follow a similar approach as in the convective problem described in Section 4.3, starting with solving the temporal instability problem to find the phase speed c corresponding to an imposed real wave number k. We then varied c away from this starting point while tracking the corresponding frequency and wave number to map out the dispersion relation (k) over a certain region in k-space. In many cases it was possible to identify a saddle point in the (complex) surface (k), as illustrated in Figure 6 for a basic state with = 6 and = 13. Note that (k) being an analytic function of the complex variable k implies that the two surfaces show saddles at the same location.
To verify that the saddle point shown in Figure 6 does indeed describe a physically relevant mode and results from the merging of modes originating from opposite sites of the real k-axis (Section 4.1), we computed the Bromwich images corresponding to a range of Bromwich contours, that is, we computed the values of the complex (inverse) dispersion relation k( ) for a range of i -isolines, as shown in Figure 7. To obtain the Bromwich contours and images we performed two sets of calculations, both initialised with a value for c close to the previously computed complex saddle-point value c a = a ∕k a , but with slightly increased or decreased (about 1%) imaginary part, respectively. For the two slightly altered values of initial c, solving Equation (3) will lead to two different solutions in k-space on either side of the saddle point (denoted as mode 1 and mode 2 in Figure 7). We then vary this initial c based on Equations (6) and (7) to obtain the Bromwich contours and the corresponding Bromwich images for the two modes individually by varying c so that the change in either i or r vanishes. Note that this procedure will generally lead to Bromwich images corresponding to different sets of Bromwich contours for the two modes, but can still provide useful information on the structure of It can be seen in Figure 7 how the Bromwich images of the two modes originate from different sides of the real k-axis for Bromwich contours with i ≫ a,i . As i → a,i for the Bromwich contours, the corresponding Bromwich images of the two modes merge at the saddle point. This behaviour confirms that the identified saddle point is indeed physically relevant and does describe a mode of absolute instability. We can perform the described instability calculation for various basic states to investigate the dependency of the absolute instability (described by the saddle point location) on the specified flow profile. Figure 8a shows the absolute frequency a and wave number k a as function of the parameter (for fixed = 6) determining the basic state according to Equation (4).
Two main observations can be made: • The temporal growth rate depends strongly on . In particular a,i changes sign at ≈ 7. Here the basic state transitions from being absolutely stable to being absolutely unstable, • the three quantities a,r , k a,r and k a,i , that is, the real parts of the absolute frequency and wave number and the imaginary part of the absolute wave number, seem not to be sensitive to changes in basic state strength. Their values only change slightly with for very weak basic flows, for which the system is absolutely stable and the shown mode is less important.
We also investigated the characteristics of the absolute instability for various values of the width of the basic state flow with fixed magnitude = 13, shown in Figure 8b. We find a similar qualitative behaviour as for changes in , in the sense that the system becomes more stable for larger length-scales and a,i does become negative for certain parameters.

NUMERICAL SIMULATIONS OF UNSTABLE GROWTH
In this section we complement the spatio-temporal stability analysis presented in the previous section by direct numerical solution of the linearised equations governing the evolution of disturbances to a specified basic state, that is, the linearised form of Equation (1).
Unstable modes of the basic state with imposed scale and magnitude are excited by perturbing the system with the mass source shown in Equation (2), with imposed perturbation frequency 0 , weak forcing magnitude F 0 = 2 ⋅ 10 −3 and radius r 0 = . In contrast to the experiments performed in Section 3, the imposed mass source is not supposed to directly force the flow, but merely to explicitly perturb the system at a given location with a given frequency. These perturbations are representing small random perturbations of the beta-plume corresponding to the steady linear solution of Equation (1). The methodological idea of the periodic forcing with fixed frequency is that it will only excite waves and unstable modes of exactly that frequency, while the small forcing magnitude ensures that only (spatially or temporally) amplifying modes lead to a significant perturbation of the system. For long times the flow perturbation should then be dominated by the corresponding fastest growing mode. Hence it should be possible to verify numerically the behaviour predicted on the basis of the dispersion relation. Note that to the far west of the forcing (x < −200) the localised sponge region is still active (although not shown in most plots), in order to simulate an infinitely extended domain. In contrast to the steadily forced experiments in Section 3, the sponge region now acts only on the zonally non-symmetric part of the response and thus does not affect the basic state. Figure 9 visualises the time evolution of the stream function response of the linearised and non-linearised version of the model to a weak time-periodic perturbation with 0 = 0.8 and basic flows with different strength parameters and fixed width = 6. This frequency 0 corresponds roughly to the absolute frequency and the frequency of maximum growth rate of the convective instability of the system (Figures 5 and 8). Specifically shown is the logarithm of the absolute value of the stream F I G U R E 9 Temporal evolution of the logarithm of the absolute value of the stream function anomaly from the basic state taken at (x = −50, y = 0) for different values of basic flow strength and using both the linear and the nonlinear model. The slope of the straight line indicates the theoretically predicted absolute temporal growth rate a,i for a system with = 6 and = 13 according to Figure 8 function anomaly from the basic state (as defined by Equation (4)) at a fixed location west of the forcing with x = −50 and y = 0.
The time period between days 0 and about 100 is associated with the spin-up of the statistically steady response. For times after day 100 the three systems show distinctly different behaviours. Of particular interest in these plots are the positions of the local maxima of the various curves, which indicate the strength of the envelope of the temporally periodic wave that is perturbing the beta-plume initial state. The case with = 6.5 shows a constant magnitude of the perturbation at the fixed location and no temporal growth, which is consistent with the idea of convective instability. The case with = 13 and the linear model, on the other hand, shows an exponential temporal growth, suggesting an absolute instability of the system. The temporal growth rate is close to the theoretical prediction a,i deduced from Figure 8 (indicated by the slope of the straight line). However, if we impose a basic state with = 13 and use the non-linearised version of Equation (1), the response shows nonlinear saturation and becomes statistically steady for long times.
Correspondingly to the temporal evolution in Figure 9, Figure 10 shows the spatial structure of the perturbations for basic states with different values for , and again using the linearised or non-linearised version of the model. Figure 10 shows two different types of plots, all corresponding to the instantaneous state of the system on day 200. The two-dimensional filled-contour plots (Figures 10c,f) show the spatial structure of the total stream function at day 200, the line curves (Figures 10a,b,d,e) show the corresponding logarithm of the absolute value of the stream function anomaly from the basic flow along y = 0.
Figures 10a,d show the response of the linear model. Both cases, the strong ( = 13) and weak ( = 6.5) basic flows, exhibit a strong growth of the spatially periodic variability in the zonal direction, indicating a spatio-temporal instability of the basic flow. In the case with = 6.5 ( Figure 10a) the spatially growing wave is (statistically) steady in a temporal sense, as was shown in Figure 9. However, the perturbations in the system with = 13 (Figure 10d) do grow in space and time and thus strongly dominate the flow after a certain time. The magnitude of the forced wave grows rapidly (in the linear case) and does not seem to be limited 1 . The persistent exponential temporal growth of the perturbations leads to extreme values of at day 200. The structure of the stream function anomaly, however, still shows spatial exponential growth, which extends continuously to the west and the east of the forcing. The temporal amplification of a spatially amplifying perturbation, and the propagation of the wave in positive and negative x-direction, both suggest an absolute 1 Note a slight deviation of the ln | − | curves (Figure 10a,d) for large negative x below the linear trend, which is likely due to the hyper-diffusion used in the model. instability of the basic flow in the case with = 13, while the temporally steady and only westward response in the case with = 6.5 suggests a convective instability of the system. The spatial growth rates of the stream function in Figure 10a,b,d are close to the predicted growth rates k i for convective and absolute instability, respectively, as deduced via Figure 5 (for = 6.5) and 8 (for = 13). The corresponding growth rates are indicated by solid green and pink straight lines. Figures 10b,c,e,f show the response of the same two basic states, but using the fully nonlinear QG model. We can observe only small differences to the linear case for = 6.5 since the perturbation of the basic flow is relatively weak and nonlinear effects are relatively small. For the state with = 13, however, the nonlinearity of the system limits the temporal and the spatial growth drastically and hence the perturbations are much smaller than when using the linearised equations. Since the instability still occurs over the full length of the domain, the initially zonally symmetric basic state now breaks into isolated vortices everywhere, in particular to the east of the forcing. In this case the system does reach a (statistical) steady and zonally almost periodic state (Figure 10f) resembling a zonal chain of discrete vortices similar to the response to a steady mass source in the shedding state 2 (e.g., Figure 1f).

F I G U R E 10
To better understand the response in terms of convective and absolute instabilities, we performed Fourier decompositions of the stream function signal (with zonal mean removed) in time and space. The computed power spectra are displayed in Figure 11. Again, there is a clear difference between the cases with = 6.5 and = 13. For = 6.5 we do find sharp temporal peaks, located almost exactly at the imposed forcing frequency 0 . Correspondingly we observe isolated and pronounced spatial Fourier peaks for the different values of 0 . Note that the Fourier peaks are normalised by their integral over the full domain range.
For the case with a stronger basic flow ( = 13), we do not observe a range of isolated temporal or spatial peaks, but both normalised power spectra are seemingly independent of the forcing frequency, showing a single robust peak in frequency and wave number. This kind of excitation 2 Recall that in the steadily forced cases, the beta-plume corresponding to the steady linear solution only extends to the west of the forcing while the experiments performed in this section are started with a zonally symmetric beta-plume initial state. Since the absolute stability analysis in Section 4.4 predicts a strong sensitivity of the temporal growth rate a,i to changes in basic flow strength and even a sign change near ≈ 7, we can expect the absolute mode to grow very slowly for states with just above this value and to not amplify at all for < 7. Figure 12 shows the temporal peaks of a Fourier decomposition over a varying time window, with fixed length of 30 days, for a constant forcing frequency 0 = 0.4 and different values of basic flow strength.
In the case with = 6.5 the system is absolutely stable, that is, a,i < 0. The response shows a clear and localised peak exactly at the forcing frequency and independent of the analysed time period. If we choose a system with = 8.6, and hence weak absolute instability, we observe the same temporal peak at = 0 when analysing early times between days 30 and 60, since the perturbation response is directly forced by the imposed time-periodic mass source. As we move the window to later times, we see the peak at the forcing frequency becoming weaker and a second peak growing at the position of the predicted absolute frequency a,r . This behaviour illustrates the domination of the system by absolute instability for later times, since the time is needed for the instability to take over scales like the absolute growth rate −1 a,i . When using a strong background flow with = 13, the response is dominated by the strong absolute instability right away, since a,i is large.
The position of the spatial and temporal Fourier peaks in Figure 11 correspond to the dominant modes of the response of the basic state when perturbed by a time-periodic forcing of frequency 0 . In the case with = 6.5, we can plot the positions of frequency peaks against the positions of the corresponding wave number peaks to get a numerically obtained dispersion relation for k r ( r ). The result is shown as triangles in Figure 5, compared to the theoretically computed convective dispersion relation of the system with = 6.5. The good match between theory and measured dispersion gives confidence that the spatio-temporal instability theory presented in Section 4 is indeed a useful explanation for the break-up of the beta plume discussed in this section, and that the developed methods can be used to theoretically predict certain aspects of the dispersion relation of the system. As argued in Section 4.4, it is not clear if the theoretical dispersion relation is valid for | r | > a,r , hence Figure 5 does not show the corresponding measured points with | 0 | > 0.8.

APPLICATION TO STEADILY FORCED EXPERIMENTS
We will now use the presented spatio-temporal instability analysis as a theory to explain and predict the different behaviours seen in the steadily forced single-layer model experiments described in Section 3.
To link the results of the theoretical stability analysis for a certain imposed flow profile and the observed response to a steady forcing, we will first define two sets of stream function profiles (which then imply e.g., the zonal wind field): • Idealised profiles: given by Equation (4) with defining parameters and . The stability analysis in Section 4 was performed for these profiles.
• Measured profiles: extracted from forced experiments (as shown in Figure 1), with steady ( 0 = 0) forcing defined in Equation (2) with magnitude F 0 and radius r 0 . The one-dimensional profiles are obtained by averaging along the x = −10 line over days 300 to 1,000. features of the measured profiles, including a reversal of the PV gradient south of the origin as well as small and smooth secondary extrema of the PV gradient near y = ± . We will perform the stability analyses described in this section only for basic states given by idealised profiles in order to simplify the subtle computation and identification of the saddle points in the dispersion relation. The zonally independent idealised and measured profiles are supposed to be different representations of the westward extending beta-plume caused by a steady localised forcing if it were not subject to any dynamic instability (e.g., Figure 1a). To link the set of measured profiles to corresponding idealised states, we fitted the stream function structures of the latter onto the former. The fit was computed by finding the parameter combination of idealised flow width and strength that minimises the square-integrated stream function difference between the idealised profile and the corresponding measured profile extracted from the experiments with steady forcing defined by F 0 and r 0 .
The idealised nature of the fit can lead to slightly changed stability characteristics compared to the actual measured profiles. Table 1 gives an example of theoretically derived stability characteristics for different measured flow profiles and their corresponding fits. As can be seen the results of the absolute stability analysis for the measured profile and the fit are fairly close to each other for the case with r 0 = 4, F 0 = 2.2(leading to a fitted idealised profile with =5.5 and =16). In particular both values of a,i predict the system to be absolutely unstable. The system with r 0 = 4, F 0 = 1.3 on the other hand shows a significant difference in stability when compared to its idealised fit (with =5.4 and =9.5). While the measured profile predicts the system to be absolutely stable, with a,i < 0, the fit suggests the opposite. Such significant differences do mainly occur for certain transition states, that is, marginally stable or unstable states with | a,i | ≪ 1, while the fitting gives reasonably good estimates of the systems stability for very stable or very unstable states.
Also note that the predicted (both, using idealised fit and measured profile) frequency of the shedding state (F 0 = 2.2, r 0 = 4) is in good agreement with the observed frequency of the steadily forced response of ≈ 0.8, as discussed in Section 3. The good agreement gives confidence that the shedding process in the steadily forced experiments can indeed be interpreted as result of an absolute instability.
Finally we can combine the gathered information about the experimental results and the results of our theoretical stability analysis of Section 5 to create a phase diagram indicating the transitions between the three presented states as we vary the width and magnitude of the forcing. For the theoretical part we expect the behaviour of the response to change in the following way: • The flow should be steady for forcing distributions that do not lead to a sign change in the PV gradient and therefore do not satisfy Equation (A3). Generally the minimum PV gradient for the idealised cos 4 profiles will become more negative (thus the system should become less stable) for stronger and narrower forcing profiles.
• For / -combinations that lead to a sign-change of the PV gradient, the flow can be unstable in a convective or absolute sense. The transition depends on the sign of the absolute temporal growth rate of the system and the parameter combination for the transition can be deduced from plots like in Figures 5 and 8. For systems with negative a,i , we expect the flow to be in the break-up state, while positive a,i (and thus absolute instability) should correspond to a flow in the shedding state.  Figure 13 illustrates the theoretical parameter combinations where either the minimum meridional PV gradient of the linear solution (solid curve) or the temporal growth rate of the absolute instability (crosses) change sign. We then plotted diagnostics derived from the steadily forced experiments on top of the lines from the idealised theory. The position of each marker corresponding to a certain set of forcing parameters r 0 and F 0 and is given by the parameters and calculated via the idealised fitting approach described earlier. The classification of the type of response for each experiment is done based on the spatial structure and temporal variability of the flow (e.g., Figure 1) as discussed in Section 3. Note that the classification will generally be difficult for the inevitable transition states discussed earlier, meaning states that are only marginally stable/unstable in a total or absolute sense.
The theoretical predictions shown in Figure 13 are in good agreement with the experimental results. Both sets of analyses show a clear division into three regions with states corresponding to break-up or convective stability being unlikely to occur for small forcing magnitudes or small length-scales. Further, the parameter combinations at which the transition between the states happens are in good agreement with the theoretical predictions from our stability analysis, keeping in mind the uncertainties associated with transition states mentioned earlier. DK89 discussed the response in experiments using a shallow-water reduced-gravity model forced by a steady localised mass source, similar to the experiments in Section 3. They restricted their analysis to a forcing with r 0 = 3 and observed only steady and shedding states, but no convectively unstable break-up behaviour. This could either be due to the fact that they used a rather short periodic channel of length L x = 60 or because the restriction of the forcing length-scale puts their experiments in a regime where the break-up state is practically impossible to observe. Figure 13 indicates three of our forced experiments with r 0 = 3 and their position within the developed phase map. As can be seen, the region with inverted PV gradient and a,i < 0 is very small here, making it unlikely to observe a corresponding break-up state. However, it is not clear what changes in convective and absolute stability are potentially caused by, for example, their different choice of̂, hence this comparison needs to be done with caution.
characteristics of unstable behaviour in either a temporal, convective or absolute sense. Unlike other authors, we did not restrict our analysis to absolutely unstable modes, but were able to also obtain information on the convective stability of the system. These algorithms can potentially be used to investigate the spatio-temporal stability of other systems. We further applied the different types of stability analysis to various sets of idealised flows and investigated their stability properties depending on different imposed flow parameters, such as the strength and meridional scale. From this analysis we argued that the steady, break-up and eddy-shedding states corresponded respectively to stability, convective instability and absolute instability of the flow forced by the mass source. Among other results, we found the studied system to develop an absolute instability for strong flows with small meridional length-scales, while weak flows with large meridional scales typically lead to convective instabilities or stable behaviour. We were able to numerically calculate various physical scales, for example, the required threshold values for the state transitions and the frequency of the resulting absolutely unstable mode for a given shear flow profile. We demonstrated good agreement between the predicted state transitions and shedding frequencies and those observed in numerical simulations.
This flow configuration studied here was originally motivated by the westward eddy-shedding phenomenon observed in relation to the Asian monsoon anticyclone, but it has also been studied as a simple model for ocean flows, such as the break-up of the Mediterranean salt tongue as it emerges into the Atlantic. Another potential application for the general theory presented in this paper might be to study the growth of waves propagating on the African Easterly Jet, which is another example of a localised thermal forcing giving rise to an unstable zonal flow. As noted previously, Diaz and Aiyyer (2015) have studied the spatio-temporal instability of this flow, focusing on the possibility of absolute instability, and it is possible that transitions between convective and absolute instability might be relevant to, for example, the seasonal evolution of the flow. While the break-up state is probably not strongly relevant to the observed behaviour of the Asian monsoon anticyclone in particular, the description of the eddy-shedding process as spatio-temporal instability serves as a new conceptual model for the phenomenon and provides quantitative predictions on scales and behaviours that have previously not been possible.

Derivation of an equivalent condition for an imposed forcing
DK89 derived an analytic criterion for which the steady linear response to a given mass source in a full shallow-water model reverses the PV gradient. We will now present the derivation of an equivalent condition for an imposed forcing F and the QG single-layer model shown in Equation (1). For simplicity, we will neglect the hyper-diffusion and the sponge region terms, that is, we set = 0 and S ≡ 0. The steady and linearised model (for a resting fluid) is then given by Similar to DK89, we can assume that only transient waves influence the region east of the forcing, giving = 0 for x → ∞, and thus solve Equation (A1) by a a corresponding integration. This leads to a formula for the meridional PV gradient: Note that y q in Equation (A2) changes monotonically with x if F is single-signed and thus its extreme values in x can be estimated for x → −∞. Under this assumption, we can also assume F ≥ 0 everywhere without loss of generality. If that is the case, a necessary condition for the meridional PV gradient to be reversed is given by Equation (A3) is equivalent to the equation DK89 found for a full shallow-water model and can be used to predict if a steadily forced system can support instabilities.