Atmospheric undular bores

We show that a recently-derived model for the propagation of nonlinear waves in the atmosphere admits undular bores as travelling-wave solutions. These solutions represent waves consisting of a damped oscillation behind a front that is preceded by a uniform breeze-type flow. The generation of such wave profiles requires a jump in the heat source across the leading front of the wave, a feature that is consistent with observations.


Introduction
The morning glory wave-structure is a quite amazing sight, but producing a robust and consistent theory for its investigation is a considerable challenge.This phenomenon comprises a travelling wave which propagates with little perturbation of the background uniform shear flow (a breeze that typically varies only with height), followed by a steep, rising front with trailing, damped oscillations.The result is a spectacular meteorological undular bore which often becomes visible during the early morning as bands of clouds, sometimes stretching from horizon to horizon (see Fig. 1).An important element in the production of these waves is the thermal inversion (the reversal of the normal tropospheric temperature reduction with height -see [6]), where a layer of cooler air is overlain by a layer of warmer air.This occurs typically after a warm airflow encounters colder, humid air.The warmer air is lifted and glides over the cold air-mass, thereby producing the thermal inversion.Then, as a wave crest approaches, the parcels of humid air rise and cool, condensation occurs and clouds are formed.Subsequently, as the wave crest passes, the parcels of air descend and warm up, so that the cloud evaporates.The clouds usually disappear by mid-morning, as the overall air temperature rises.It should be emphasised that, in this process, the clouds are not carried along with the air flow; rather, clouds are continuously being formed at the leading edge of each wave and continuously eroded at the trailing edge.Atmospheric undular bores often propagate at speeds in excess of 10 m/s and may last for hours (see [11]).More details of the background to this phenomenon, and to the underlying physics, can be found in [8].
As we have indicated, the mathematical description of these atmospheric undular bores is far from routine work.It is a fairly straightforward exercise to write down model equations that appear to possess the properties observed in the morning glory; see, for example, [4,12,18].Some of the attempts at modelling involve arguing that a direct connection exists between these waves and the Korteweg-de Vries or Benjamin-Ono equations; see [2].Typically, the time and length scales in these model equations are not quite correct but, by analogy, it is argued that something along these lines is the appropriate starting point.The application to the morning glory waves then requires the scales and available coefficients to be matched to the data in order to recover what is observed; see, for example, [14,22].Of course, numerical simulations provide some useful clues, as does the application of arguments that draw an analogy with the undular bores observed in rivers.We refer to [15,16] for a discussion of hydraulic bores but we note, even in this simpler context, that it is difficult to construct closed-form representations of the wave profiles.Moreover, whereas hydraulic bores in water are typically generated by fast tidal currents being funnelled into narrow estuaries (see [1] for some beautiful photographs of tidal bores on the Severn river), atmospheric bores are formed by warm air flowing over colder air.Also, in stark contrast to water flows, which typically transport energy but not mass (see the discussion in [5]), atmospheric bores transport air and energy (see [19]).For example, the lifting of the air parcels by the passage of the bore, with nett vertical displacements of about 1 km, plays an important rôle in the nocturnal convection over the Great Plains in the central US (see [11]), which produces most of the warm-season precipitation [17].While high humidity makes the magnificent morning glory clouds over the Gulf of Carpentaria in Australia clearly visible (see [20]), bores over the central US are typically traced by means of radar images (partly because they often occur during the night, when direct confirmation by an observer is difficult and satellite imagery is not available).
In all the material cited above, and in the many others recorded in [8], there is no careful and comprehensive derivation of a suitable equation (or system) from the governing equations that describe the fluid dynamics of the atmosphere.That this is possible is demonstrated in detail in [8].(We provide a brief overview of the relevant derivation, so that the reader is able to put the results of the current work directly in context.)A number of special solutions were highlighted in this earlier paper, each extracting an important property of these wave systems, but very little more was offered.The work presented here attempts to combine all the salient features of the morning glory in one solution of this model equation.As we shall see, this new solution captures all the properties of atmospheric undular bores.The plan for this paper is as follows.After presenting in Sect. 2 a brief outline of the derivation of the equation which describes nonlinear wave propagation in the troposphere, we show in Sect. 3 how a Liénard-type second-order ordinary differential equation underlies the modelling of undular bores.In Sect. 4 we investigate the qualitative behaviour of the solutions of the Liénard equation; this provides a firm basis on which to present a discussion of the wave structure that is typical of the undular bores.Finally, in Sect.5, we combine all these results to present a detailed mathematical description of the morning glory, which includes the rôle of the background atmosphere and associated breezes, as well as a numerical simulation that depicts a typical undular bore.

Preliminaries
We rely on a recently-derived model for the nonlinear propagation of atmospheric waves (see [8]), derived from the governing equations for flows in the troposphere, written in rotating spherical coordinates -the Navier-Stokes equation with variable eddy viscosities, the equation of mass conservation for a compressible fluid, the equation of state for the atmosphere and the first law of thermodynamics -using only the thin-shell approximation.The model accommodates waves propagating horizontally through the atmosphere in any direction, with the underlying asymptotic structure requiring that the length scale of the wave packet should be the same as the thickness scale of the troposphere (about 16 km), with the length scale of the wave front based on the radius of the Earth.Thus the wave packet might extend tens of kilometres, but the wave front could be a thousand kilometres long.The time scale of the motion was given by 1/Ω ≈ 3 1  2 hours, where Ω is the constant angular speed of rotation of the Earth.
The spherical coordinates (ϕ, θ, r ) in the rotating frame are the distance r from the centre, the angle of latitude θ ∈ [− π 2 , π 2 ] and the angle of longitude ϕ ∈ [0, 2π).They are redefined as for a clockwise rotation of the coordinates about the fixed point (ϕ 0 , θ 0 ) through the fixed angle α; here ε = H /R is the thin-shell parameter, with H the maximum thickness of the troposphere and R an average radius of the Earth.Throughout, the prime denotes dimensional (physical) variables so, in particular, we write the time variable as t = t/Ω and the radial variable as r = R (1 + ε Z ), where t and Z are nondimensional variables.Further, the horizontal velocity components, (u , v ), in this spherical system, are rotated to produce components that are in the direction of propagation (V ) and along the wave front (U ): the vertical component, w , is unchanged; we use H Ω ≈ 1.2 m/s as the speed scale.Further details of the non-dimensionalisation are given in [8], and hereafter nondimensional variables are written without the prime.The process then involves the perturbation of a general background state of the atmosphere, eventually producing a nonlinear equation for the stream function, , which describes the leading-order approximation, as ε → 0, of the velocity components (V 0 , w 0 ) with the leadingorder component along the wave front taken as U 0 = 0.The evolution equation for ( , , Z , t), at leading order, is where ρ 0 (ξ ) with ξ = gZ + 1 2 S 2 is the background density in the atmosphere, and with g ≈ 9.8 m/s 2 the average acceleration of gravity at the surface of the Earth.We have also introduced: an average density of the atmosphere ρ and an average dynamic eddy viscosity μ , with μ μ v (Z ) and μ μ h (Z ) the vertical and horizontal dynamic eddy viscosities, respectively.The Reynolds number, R e , has been treated as O(1) under the limiting process ε → 0, so that the dominant viscous effects are retained.The forcing term, K ( , Z , t), is any appropriate external heat-forcing that drives and maintains the motion.We note, in passing, that any variation along the wave front, which depends on , appears only parametrically in this equation.Finally, we record the velocity components expressed in terms of the stream function: Of special interest are travelling-wave solutions that describe wave propagation for fixed , at constant speed c( ), without change of shape.These solutions correspond to stream functions of the form and associated velocity components In the two-dimensional frame of reference moving horizontally at the wave speed c( ), the particle paths t → ( , (t), Z (t)), obtained by solving the system (since U 0 = 0 ensures that remains constant for each individual particle), coincide with the level sets of the function because (2.3) and (2.4) yield This shows that the profile of the stream function , at fixed , offers insight into the global dynamics of these atmospheric flows.A simplification of the explicit appearance of the density ρ 0 in Eq. (2.1) can be obtained by introducing This is valid for any general form of the density throughout the troposphere, and the positivity of the density ensures that we may treat ρ 0 = d(y, z), in terms of the new variables (x, y, z) with ( , ) = (x, y).Setting = ψ(x, y, z, t), Eq. (2.1) becomes (2.6) where Correspondingly, the velocity components (2.2) can be written as Three important types of solution to Eq. (2.6) were obtained in [8], for special choices of the heat forcing k, and of the vertical and horizontal eddy viscosities m and M: • Breeze-like horizontal flows (uniform in the x-variable at every fixed height), described by stream functions of the form ψ(z, t).Such flows are typically initiated by the different solar heating of land and water surfaces, with marine air advancing onto land as a sea breeze (see [23]).• Shallow-water bore-like solutions, described by stream functions of the form where χ > 0, κ > 0, c > 0 depend only on y, representing smooth travellingwave fronts that bring about the transition between the levels b(y, z){1 + χ } for x → −∞ and b(y, z){1 − χ } for x → ∞ (relative to a).This simple model of a bore is commonly used in shallow-water theory -a framework unable to predict undulations because the vertical component of the velocity is neglected (see the discussion in [9]).Note that the vertical velocity component in an atmospheric undular bore often exceeds 2 m/s (see the data in [11]).
• Uniformly oscillating solutions, described by stream functions of the form ψ(x, y, z, t) = a(y, z) + χ e κz cos(κ[x − ct]) where χ = 0, κ = 0, c = 0 depend only on y.Such motions are typically observed as adiabatic buoyancy oscillations in a stably stratified atmosphere and are also captured by linear theory (see [7]): for average tropospheric conditions the period of such buoyancy oscillations is about 8 min [13].However, damped oscillations are a hallmark of non-turbulent atmospheric bores (see the discussion in [24]).
In this paper we present a solution which essentially combines all these elements (a breeze-like flow as a precursor of a bore, with the jump followed by damped oscillations) and so confirms the relevance of Eq. (2.1) to the modelling of undularbore phenomena in the atmosphere.Our approach goes beyond the usual analyses that use, piecemeal, different models to describe different components of the atmospheric undular bore.We show that this new travelling-wave solution captures all the salient features of phenomena such as the morning glory: damped oscillations behind the leading front, with a horizontally uniform flow ahead and a wave structure generated by a jump in the heat source across the front -a feature consistent with observations.

Travelling-wave solutions with damped oscillations
The over-arching assumption that we make here is that the solution of interest is expressed as a height-dependent travelling wave in the form where c is the constant speed of propagation; the parametric dependence on y has been included, but it will play no rôle in the development, so we will suppress it hereafter.(It could be reinstated and used, for example, to limit the extent of the wave front in the y-direction.)For (3.1) to be a solution of (2.6) requires where p, q and r are constants.These three equations are in the two unknowns, a and b; we must therefore accept that we shall have to make special choices for m and/or M in order to complete the solution.We anticipate, however, provided that we retain some appropriate form of the variable viscosities, that the fine detail of how the viscosities vary will not significantly alter the main structure of the solution.The resulting equation for F(x − ct) now becomes where ξ = x − ct and k 0 is a constant.Note that: • Setting p = 0 in (3.5) leads to no simplification, since the new form of (3.5) can be brought back to the original version by adding suitable constants to F and k 0 .• For q = r = 0 Eq.(3.5) simplifies to the separable equation • For r = 0 and q = 0 Eq.(3.5) simplifies to with general solution given implicitly by for some ξ ∈ R.
• For q = 0 and r = 0 Eq.(3.5) simplifies to transformed by integration to the Riccati equation for some k ∈ R. The substitution F = −2r ϕ /ϕ transforms it into the linear second-order differential equation For k 0 = 0 the general solution of the corresponding constant-coefficient equation is for arbitrary constants ϕ 1 , ϕ 2 ∈ R. For k 0 = 0 the general solution is given by (see [21]) where J 1/3 and Y 1/3 are the Bessel functions of the first and second kind, while ϕ 1 , ϕ 2 ∈ R are arbitrary constants.Since for a non-trivial solution ϕ of (3.6) the functions ϕ and ϕ cannot vanish simultaneously, we see that the above oscillatory solutions to (3.6) correspond to solutions of (3.5) that blow-up.
The previous discussion shows that for q = 0 and r = 0 equation (3.5) does not admit oscillatory solutions.
For qr = 0 it is convenient to normalise Eq. (3.5) by writing where F 0 = 0 and ξ 0 > 0 are constants; we then choose Under these transformations, the velocity components become Before we embark on the construction of a solution of Eqs.(3.2)-(3.4)and (3.7), we make a few general observations about Eq. (3.7), which is the equation at the heart of the analysis.This equation looks disarmingly simple, but it is not!To see this, we transform according to which gives the Abel equation of the second kind where Y = 1 2 ( f + ν) 2 and the choice of ± in (3.9) is given by the sign of ( f + ν).The fact that the forcing term in (3.9) combines a constant with a term proportional to 1/ √ Y immediately precludes any possibility of finding a solution in closed form, unless κ + ν sgn(rq) = 0 (see the discussion in [21]).For κ + ν sgn(rq) = 0 we can proceed, as (3.9) simplifies to but this is not a relevant choice for solutions of interest here -there are no oscillatory solutions.Despite the non-availability of explicit solutions for κ + ν sgn(rq) = 0, we are able to make a simple observation which allows us to make some headway in this case.We seek a solution f to (3.7) which approaches a constant f 0 as X → −∞ and which oscillates (this being an appropriate behaviour for the upper level in an undular bore); to do this we set where f 0 , s and λ are constants, with Re(λ) > 0 and Im(λ) = 0.This is a suitable asymptotic solution of Eq. (3.7) when we have and then the solution that we seek requires This asymptotic result can be embedded within a more rigorous analysis of the dynamics of the equation as we show in the next section.Eq. (3.11) is obtained by taking sgn(rq) = 1 (as we found above) and setting

Global dynamics of the Liénard equation
We now discuss the global dynamics of the autonomous Liénard Eq. (3.11) with parameter η ∈ (1, 2).There is a vast research literature devoted to qualitative studies of Liénard type equations (see the discussion in [25]) but mainly from the perspective of limit cycles; the specific form (3.11) is not encountered in these considerations.Equation (3.11) looks disarmingly simple, but, as discussed in Sect.3, there are no explicit non-trivial solutions available.We prove the following result.
strictly increasing as k ∈ Z increases.
we see that Theorem 1 ensures that a specific solution of (4.1) spirals in towards the unique equilibrium point at the origin for X → −∞, and spirals out from it for X → ∞; see Fig. 2. Since by uniqueness the solution curves cannot cross, this, in combination with the behaviour of the vector field along the two nullclines (see below), shows that every single non-trivial solution curve of (4.1) traces a similar spiralling pattern.
Proof To establish the solution behaviour claimed in Theorem 1, it is convenient to study (3.11) for X < 0 by reversing the flow direction.More precisely, setting we transform (3.11) for X < 0 into Since the characteristic equation for the linearisation of the system (4.1)near the origin is , and that of the linearisation of the system (4.4) near the origin is , the Hartman-Grobman theorem (see [3,10]) applies.Thus, since 1 < η < 2, the origin is a repelling focus for the nonlinear system (4.1) and an attractive focus for the nonlinear system (4.4).This confirms, in part, the behaviour of the solution claimed by Theorem 1, but the detailed global pattern specified therein requires a more detailed investigation, based on an interplay of phase-plane analysis and the construction of suitable Lyapunov functions.It turns out that the analysis for (4.4) is rather routine, whereas tackling (4.1) is quite intricate.
The vertical nullcline of (4.4) is the parabola [g = η f − 1 2 f 2 ], along which the vector field points vertically upwards in the left half-plane [ f < 0] and vertically downwards in the right half-plane [ f > 0].The horizontal nullcline of (4.4) is the axis [ f = 0], with the vector field pointing to the right along the positive semi-axis and to the left on the negative semi-axis.Since the origin is the only equilibrium point of (4.4), and solution trajectories are disjoint by uniqueness, this information suffices to infer that the claimed spiralling behaviour of the nontrivial solutions to (3.11) for X < 0 follows once we establish it for a single solution.With this aim, we now investigate the behaviour of the unique solution curve C(X ) = ( f (X ), g(X )) of (4.4) with C(0) = (2η, 0), that is, with f (0) = 2η and g(0) = 0 (4.7) so d f dX (0) = 0 , defined on its maximal existence interval (X * − , X * + ), with X * + > 0 and X * − < 0 possibly finite; see Fig. 3. Using the Lyapunov function from (4.7) we deduce that the interior O of the circle [ . Indeed, the solution curve C enters this region since (3.11) and (4.4) yield 2 for all ε > 0 small enough, due to the Taylor expansion Also, once inside the region O, throughout which f 2 < 4η 2 (and thus −2η < f < 2η), the relation (4.9) ensures that X → L( f (X ), g(X )) decreases with growing X .In particular, the solution remains bounded and thus blow-up at finite X > 0 is not possible, so X * + = ∞.Since in the region [ f > 0, g < 0] ∩ O, where by (4.10) the solution curve X → C(X ) is located for X > 0 small enough, the functions X → f (X ) and X → g(X ) are strictly decreasing, due to (4.4), with X → f 2 (X ) + g 2 (X ) decreasing by (4.9), we deduce that the solution curve hits first the negative semi-axis [ f = 0, g < 0] at some X + 0 > 0, with Fig. 3 The solution curve X → C(X ) of the nonlinear system (4.4) through (0, 2η) spirals around the origin (the single equilibrium point of the system), with continuously diminishing distance from it as X grows towards ∞.The vector field points vertically upwards along the dotted branch of the vertical isocline in the half-plane [ f < 0], vertically downward along the dashed branch in the half-plane [ f > 0], and horizontally along the vertical axis f = 0 (to the right on the dotted positive semi-axis and to the left on the dashed negative semi-axis).The interior O of the disk [ f 2 + g 2 = 4η 2 ] is positively invariant.The connection to Fig. 2 is made by direction reversal, with X 1 and X 2 corresponding to X −1 and X −2 in Fig. 2 Now (4.4) yields d f dX (X + 0 ) < 0 and dg dX (X + 0 ) = 0, so the solution curve will enter the region [ f < 0, g < 0] ∩ O for X > X + 0 small enough.It is impossible for the solution curve to remain below the parabola [g = η f − 1 2 f 2 ] for all X > X + 0 since then X → f (X ) would be strictly decreasing on [X + 0 , ∞) and would make g(X ) positive for X > X + 0 + 1 large enough.Thus we can find X 1 > X + 0 so that the solution curve intersects [g = η f − 1 2 f 2 ] with f (X 1 ) < 0 and g(x 1 ) < 0, crossing the parabola vertically upwards as d f dX (X 1 ) = 0 and dg dX (X 1 ) = − f (X 1 ) > 0. The solution curve therefore enters the region Throughout this region, (4.4) yields that X → f (X ) and X → g(X ) are strictly increasing.The solution curve cannot reach the origin at finite X (since then, by uniqueness, it would have to be stationary for all values of X ∈ R) and it also cannot reach the origin as X → ∞ while confined to this region (since the Hartman-Grobman theorem ensures damped oscillations of f (X ) close enough to the origin).Consequently there exists X + 1 > X 1 where the solution hits the negative semi-axis: g(X + 1 ) = 0 and 0 > f (X + 1 ) > g(X + 0 ), with the inequality on the right side ensured by the fact that X → L( f (X ), g(X )) is strictly decreasing on (X + 0 , X + 1 ).From (4.4) we see that d f dX (X + 1 ) > 0 and dg dX (X + 1 ) > 0, so that the solution curve enters the region [ f < 0, g > 0] ∩ O for X > X + 1 small enough.In this region, (4.4) ensures that X → f (X ) and X → g(X ) are strictly increasing.Again, the origin cannot be reached in this region at finite X or for X → ∞, so that there is some 1 where the solution curve hits the positive vertical semi-axis: f (X ++ 1 ) = 0 and 0 < g(X ++ 1 ) < − f (X + 1 ), with the last inequality a consequence of the monotonicity of X → L( f (X ), g(X )) on (X + 1 , X ++ 1 ).From (4.4) we now see that d f dX (X ++ 1 ) > 0 and dg dX (X ++ 1 ) = 0, so that the solution curve enters the region 1 small enough.Throughout this region the solution curve must intersect the vertical nullcline at finite X since otherwise (4.4) would yield that X → f (X ) is strictly increasing and X → g(X ) is strictly decreasing on (X ++ 1 , ∞), but the origin is the only possible accumulation point in the closure of O. Let X 2 > X ++ 1 correspond to the intersection of the solution curve with the vertical isocline.Then d f dX (X 2 ) = 0 and dg dX (X 2 ) < 0, so that the solution enters the region for X > X 2 small enough.In this region, (4.4) ensures that X → f (X ) and X → g(X ) are strictly decreasing, and since, as above, the solution curve cannot reach the origin in this region (at finite X or for X → ∞), there is some X + 2 > X 2 where the solution curve hits the positive horizontal semi-axis: g(X + 2 ) = 0 and 0 < f (X + 2 ) < g(X ++ 1 ).From (4.4) we now get d f dX (X + 2 ) < 0 and dg dX (X + 2 ) < 0, so the solution crosses the positive horizontal semi-axis, entering the region [ f > 0, g < 0] ∩ O for X > X + 2 small enough.We now repeat the arguments step-by-step, showing that the solution spirals around the origin, forever getting closer to it as X increases.Thus we have proved the claim regarding the behaviour for X < 0 of the nontrivial solution to (3.11) with initial data corresponding to the initial data (4.7) for (4.4).
Let us now investigate the behaviour on the maximal existence interval [0, −X * − ) of the solution ( f (X ), g(X )) to the planar system (4.1) with initial data f (0) = 2η and g(0) = 0 , (4.12) corresponding to the initial data for (3.11).Rather than working with the system (4.1), it is convenient to transform it by means of f Note that (4.15) and (4.16) yield These considerations, in combination with (4.18), show that the solution curve starting at the point (−2η, 0) on the boundary of the disk [ f 2 + g 2 = 4η 2 ] is confined to the exterior of this disk for all X ∈ (0, −X * − ).From (4. 19) we see that the solution curve enters the quadrant [ f < 0, g > 0] for X > 0 small enough, and throughout this quadrant X → f (X ) and X → g(X ) are both strictly increasing, with the second statement a consequence of the fact that, in this region, any point ( f , g) with f 2 + g 2 > 4η 2 lies above the graph of the parabola g = − 1 2 f 2 − η f (see Fig. 4).We claim that there is some X † 0 ∈ (0, −X * − ) where the solution curve hits the positive vertical semi-axis: Indeed, otherwise, by choosing ε > 0 small enough, we could ensure that But then integration of (4.18) on [ε, X ] for X ∈ (ε, −X * − ) prevents blow-up at finite X , so −X * − = ∞, which is incompatible with the solution curve remaining in this quadrant since then so that integrating the first equation in (4.15) on [ε, X ] would ensure that f (X * ) = 0 for some X * > ε.
At X = X † 0 we get from (4.15) that d f dX (X † 0 ) > 0 and dg dX (X † 0 ) = 0, so that the solution curve enters the quarter-plane [ f > 0, g > 0] for X > X † 0 , with its distance from the origin increasing with growing X , due to (4.18).In this quarterplane X → f (X ) is increasing while X → g(X ) decreases.We claim that there is some X † † 0 > X † 0 where the solution hits the positive horizontal semi-axis: when we also take (4.18) into account.Indeed, to see that blow-up at finite X is not possible in this quarter-plane, we consider the Lyapunov function

.22)
If the solution curve were to remain in the quarter-plane [ f > 0, g > 0] for all X > X † 0 , then 0 < g(X ) ≤ g(X † 0 ) for all X > X † 0 and (4.22) would give so that the Gronwall inequality would yield the boundedness of L 2 ( f (X ), g(X )) on the finite interval [X † 0 , −X * − ).Knowing already that g(X ) is bounded, f (X ) would also be bounded, thus preventing blow-up.We can also rule out the possibility that the solution is defined for all X > X † 0 and the solution curve remains in the quarterplane [ f > 0, g > 0] for all X > X † 0 , since the monotonicity of X → f (X ) ensures f (X ) ≥ f (X † 0 + 1) > 0 for all X > X † 0 + 1, so that by integrating the second equation in (4.15) we would obtain Fig. 4 Behaviour of the solution to (4.15) with initial data (4.16).The connection to Fig. 2 is made by means of the transformation (4.14) and noticing that X † 1 and X † 2 correspond to X 1 and X 2 in Fig. 2, respectively and this means that g(X ) must vanish at finite X .Consequently there is some where the relations specified in (4.20) hold.From (4.15) we now see that d f dX (X † † 0 ) > 0 and dg dX (X † † 0 ) < 0, so that the solution curve enters the quarter-plane [ f > 0, g < 0] for X > X † † 0 , getting, by (4.18), further away from the origin with growing X .We now show that the solution cannot blow-up at finite X while remaining in this quarter-plane.Indeed, since the function g → g e g is bounded on (−∞, 0], Gronwall's inequality ensures, by means of relation ( 4 /2 is also bounded, X → f (X ) has to be bounded on a finite interval [X † † 0 , −X * − ), preventing the assumed blow-up at finite X .
We now claim that the solution curve must intersect the vertical nullcline We now introduce the Lyapunov function if (4.23) holds, since then the first equation in (4.15) would ensure, in view of the assumption (4.23), that X → f (X ) is increasing for X > X † † 0 , and f (X and the first equation in (4.15) would lead us to Dividing both sides of the above differential inequality by f 2 (X ) and integrating, we would obtain the inequality which clearly fails for X > X † † 0 large enough.Consequently there is some X From (4.15) we see that in the quarter-plane [ f > 0, g < 0] the solution curve crosses the vertical isocline vertically downwards, so that it moves, for in which both X → f (X ) and X → g(X ) are strictly decreasing.For ε > 0 small enough we then have with the monotonicity properties ensuring that g(X ) 1 + ε for which the solution lies in the region R. From (4.15) we now infer d f dX (X ) ≤ −E and 0 > dg dX (X ) > − f (X † 1 + ε) for all X > X † 1 + ε for which the solution lies in the region R, so that blow-up at finite X is prevented and the solution curve must intersect the negative vertical semi-axis at finite X : there exists some Since (4.15) yields d f dX (X † † 1 ) < 0 and dg dX (X † † 1 ) = 0, the solution curve enters the region throughout which, due to (4.18), it satisfies , with X → f (X ) strictly decreasing and X → g(X ) strictly increasing.Consequently the solution curve must hit the vertical nullcline at some 2 ) = 0 and dg dX (X † 2 ) > 0, so that the solution curve crosses the nullcline vertically upwards at X = X † 2 .Since the vector field points vertically upwards all along the vertical nullcline in the quarter-plane [ f < 0, g < 0], we see that the solution must lie above the curve [g + 1 2 f 2 + η f = 0] for X > X † 2 , as long as f (X ) < 0 and g(X ) < 0. Therefore, the solution curve must hit the negative horizontal semi-axis at some X † † 2 > X † 2 , with g(X † † 2 ) = 0 and f (X † † 2 ) < g(X † † 1 ) < −2η.
From (4.15) we see that d f dX (X † † 2 ) > 0 and dg dX (X † † 2 ) > 0, so that the solution curve enters the quarter-plane [ f < 0, g > 0], where it cannot cross the solution path ( f (X ), g(X )) 0<X <X † 0 , so that in this region we have d f dX (X ) > 0 and dg dX (X ) > 0 along the solution curve.We can now repeat the arguments used before to infer that the solution curve hits the positive vertical semi-axis, then enters the quarter-plane [ f > 0, g > 0] and goes around the origin, exactly as for the case X ∈ [0, X † † 2 ], but now further away from the origin, due to (4.18).We conclude that the solution curve exists for all X > 0, and spirals around (and away) from the origin, as claimed.
Thus we have a decaying oscillation in X > 0 and, with the same initial conditions for Eq.(3.11), a growing oscillation in X < 0; this result must now be embedded within a description that corresponds to the observed properties of the morning glory.It is immediately clear that we cannot permit the growing solution, but we could treat the decaying oscillation as starting from any trough (at a finite point) and evolving to the left; the solution to the right must take a different form.How this is accomplished we describe in the next section.

The undular bore
The first stage in the construction of a solution is quite straightforward: we make choices for a(z) and b(z) in (3.due to (3.5).We now turn to the all-important issue of finding suitable functions F for stream functions of type (3.1).
The considerations in Sects.3 and 4 make clear that the construction of a suitable function F will require some care.Indeed, a description of the whole flow field is a significant challenge, but showing that an undular-bore profile exists is attainable.To do this we concentrate our efforts on the streamline that defines the front of the wave and its associated oscillations.Let the streamline that is this profile be given by ψ(ξ, y) = 0, then the equation for its shape is written in terms of z and X = ξ/ξ 0 = (x − ct)/ξ 0 , where F 0 = 0 and ξ 0 > 0 are constants (introduced in Sect.3).The function which controls the shape of the wave, f (X ), satisfies see Eq. (3.7), where we have set sgn(rq) = 1.We seek a solution of Eq. (5.6) for which f → 0 as X → ∞ (ahead of the wave front) and f → f 0 > 0 as X → −∞ (behind the front, where oscillations die out).Let the corresponding solutions of (5.5) be z = z ± , so that Further, consider a point (X , z) = (X 0 , z 0 ) on the leading face of the bore, and then see Fig. 5.The analysis performed in Sect. 4 shows that a solution to (5.6) which oscillates about f = f 0 > 0 for X → −∞, that is, with κ = f 0 , cannot satisfy decay conditions ahead of the front (as X → ∞); conversely, a solution which does decay ahead (so κ = 0) cannot exhibit decaying oscillations about f = f 0 > 0 as X → −∞.We must, perforce, construct a solution which satisfies different versions of (5.6) in different regions.
In the light of the situation described above, we consider two regions: X > X 0 and X < X 0 ; we expect some form of discontinuity across X = X 0 but, as we shall see, this can be regarded as consistent with the general framework that corresponds to the generation of morning-glory waves.Thus we examine, for f 0 > 0, the equation The results in Sect. 4 show that in the region X < X 0 there are solutions (which depend on the parameter ν, but with −2 < f 0 + ν < 0) which satisfy at some X ∈ R, and which oscillate around f = f 0 and decay as X → −∞.
Correspondingly, there is the solution f ≡ 0 in X > X 0 ; it is convenient hereafter to set X 0 = 0 (see Fig. 5).Thus the profile described by (5.9), with the conditions just outlined, is a continuously differentiable function on R. The jump discontinuity of d 2 f dX 2 across X = X 0 = 0, which relates directly to the heat-source term, is consistent with what is observed as the morning glory is generated.The front of the wave is at the head of a warm breeze that moves into a colder region of air, riding over a cold breeze: there is a difference in temperature and heating across the front.In our solution, this arises by virtue of the discontinuity in d 2 f dX 2 because this is given by the jump in the forcing constant (from 0 to f 0 ).In terms of the associated stream function (3.1), this jump discontinuity of d 2 f dX 2 across X = X 0 = 0 corresponds to a Dirac-delta singularity of the vorticity, concentrated on a vortex patch.Now that we have a solution for f (X ) which represents an undular bore, and so recovers the essential structure of the morning glory, we can examine the nature of the velocity field associated with this wave.Although many different avenues and choices could be explored, we consider the simple case described by (5.1) with associated constant values of M > 0 and m > 0 in (5.2)-(5.3),while the velocity components are given by (2.7).We impose w 0 = 0 on z = 0 (i.e., no flow through the surface of the Earth, so z = 0 is a streamline): b(0) = 0 and so b (5.10) and we see that w 0 = 0 whenever d f dX = 0, which occurs in X > 0 and for specific values of X < 0; otherwise, w 0 oscillates as d f dX oscillates.The horizontal velocity component along z = 0 can also be found: A breeze exists along z = 0 because our choices do not admit a no-slip condition; let it be given by V 0 (0, X ) = V + for X > 0, V − for X → −∞.where, from our earlier considerations, −2 < f 0 + ν < 0. Thus (5.12) The colliding breezes are defined by the initial state of the atmosphere, before the formation of the morning glory; thus we may presume that we are given V + and V − .The constants C, α and f 0 are known from the prescribed state of the atmosphere, and so (5.12) determines the speed of the wave, in terms of the parameter ν ∈ (−2− f 0 , − f 0 ); an example of an undular-bore solution of Eq. (5.9) is shown in Fig. 6.The importance of the free parameter ν is revealed by interpreting (5.12) as a relation between the wave speed and the oscillation amplitude of the bore: knowing c, we obtain ν and the unique solution of the ordinary differential Eq. (5.9) in X < 0, with initial data f (0) = d f dX (0) = 0, which determines the oscillations, the size f 0 of the jump being prescribed by the height of the warmer-breeze layer on top of the colder layer which is at ground level.

Discussion
A consistent theoretical model for atmospheric undular bores was only recently developed (see [8]), while earlier studies rely on rather ad hoc simplifying assumptions, reasoning by analogy with shallow-water theory -an approach that fails to address the inherent complexity of the underlying thermodynamical processes.Another shortcoming of earlier attempts to understand these waves is the use of piecemeal models to describe, separately, the core features of the flows (such as jumps, oscillations and the existence of background breezes).So, although these methods have led to the identification of the main ingredients that make-up the morning glory, complete and relevant solution has been found that combines all the essential elements.In the current work, we remedy this situation and show that a suitable solution does exist.The travelling-wave Ansatz leads us to a nonlinear second-order ordinary differential equation -of Liénard type -that looks deceptively simple but has no explicit non-trivial solutions.By means of a quite intricate global analysis in phase-space, which relies on three functionally independent Lyapunov functions, we show that all nontrivial solutions are oscillatory, bounded (and approaching a constant state) in one direction, but unbounded in the other.This confirms the main structure of a damped oscillation at an upper level -an essential ingredient for an atmospheric undular bore -but the streamline which represents the front of the wave should also feature non-oscillatory decay ahead.Consequently different versions of the Liénard equation are required ahead and at-and-behind the front.The discussion of the global dynamics has proved that the oscillating solution exists behind and, emanating from any chosen trough, we may have a constant solution ahead.However, although the resulting solution is continuously differentiable, it possesses a jump in the second derivative at the junction of the two solutions.In the context of the morning glory such a discontinuity is consistent with the underlying physics: the jump in the second derivative comes about by virtue of a jump in the heat-source term.This describes precisely what is observed as a morning-glory wave is generated, that is, a layer of warm air over-rides a layer of colder air and so, at the front, there is a jump in temperature (and in the associated heating).This work, we suggest, therefore lays the foundations for more intensive mathematical investigations as well as careful applications to, and interpretation of, observational data.

Fig. 2
Fig. 2 Sketch of the oscillating behaviour of the solution X → f (X ) in Theorem 1.The figure on the left corresponds to the spiralling pattern of the solution curve X → C(X ) = ( f (X ), g(X )) depicted on the right in the phase-space

. 15 )
with initial data f (0) = −2η and g(0) = 0 .(4.16) Analogous to (4.8), a useful Lyapunov function for (4.15) is 1), using (3.2)-(3.4).The simplest way to proceed is to set b(z) = b 0 + b 1 e −Bz , (5.1) where b 0 , b 1 and B are constants (and we require B > 0 since at great heights -near the tropopause -this type of atmospheric flow is barely noticeable); other choices are possible, but this one is computationally advantageous.Equation (3.4) then gives M(z) = rb 0 B , (5.2) which specifies the horizontal kinematic eddy viscosity.Equation (3.3) now produces an expression for the vertical dynamic eddy viscosity m(z) = A e Bz + (qb 0 − B)b 1 B 2 , (5.3) where A is an arbitrary real constant.The various constants here can be chosen to model the viscosity, at least over the vertical extent of the undular bore; in particular, the choice A = 0 gives m constant.Finally, we use equation (3.2) to obtain an expression for a(z), to produce the associated heating source in the form k(z) = −k 0 b 0 b 1 e −Bz + β

Fig. 5
Fig. 5 Sketch of an undular bore, generated when a thermal inversion occurs -with the layer of colder air (depicted in blue) overlain by a layer of warmer air (depicted in orange) = c + pb 0 B, db dz (0) = b 0 B, due to (5.10) and (3.2).