Finding the Fibonacci Spiral Patterns in Atmospheric Equatorial Waves

Ranis Ibragimov, Eleanor Blake

Wenatchee Valley College, Wenatchee, WA, USA.

DOI: 10.4236/jamp.2022.1012250                

Abstract

The purpose of this paper is to analyze and visualize the exact invariant solution of the nonlinear simplified version of the shallow water equations which are being used to simulate equatorial atmospheric waves of planetary scales. The method of obtaining the exact solution is based on the Lie group invariance principle. It is shown that the obtained invariant solution has a Fibonacci spiral-like form and has two parameters k and t₀. We have defined a new model hypermarameter Δ_kt = t – t₀, where t is time. The question of particular interest is: can we tune the hypermarameter in order to match the exact solution to the actual Fibonacci spiral? It was discovered that the physically relevant part of the solution matches exactly the Fibonacci spiral.

Keywords

Atmospheric Waves, Fibonacci Sequence, Shallow Water Theory, Exact Solutions

1. Introduction

The Fibonacci sequence of numbers $F_n$ is defined using the recursive relation $F_n=F_{n-1}+F_{n-2}$ with the seed values $F_0=F_1=1$. It is known as probably the most common relationship that occurs in the smallest, to the largest objects in nature. It can be seen in the arrangement of sunflower seeds, shells of snails, in the pattern found on the wings of dragonflies, in the arrangement of hexagons in the honeycomb of bees, as well as in the shape of the Milky Way or other galaxies (see Figure 1). In biology, for example, the packaging of seeds or branching is only a glimpse of underlying heavy usage of Fibonacci sequence in whole morphogenetic process, starting with DNA itself. It is used in shape itself when necessary but not used when not, e.g. our limbs don’t grow in Fibonacci rule way. But it’s hidden somewhere in the shape of our body, bones etc. As a part of

Figure 1. Fibonacci spirals observed in nature. These images were posted in Art/Science, Paul’s Journal, Workshops. Available at: https://tumamocsketchbook.com/category/journal.

chemical activity, segmented spiral waves were also found in reaction-diffusion systems ( [1] [2]). The segmented spirals are complex dynamical structures with several levels of organization, characterized, for example by the wavelength between successive spiral turns and the turning wavelength between segments. The appearance and complexity of segmented spirals are suggestive of phenomena in living systems.

While the aesthetics and symmetry of Fibonacci spiral patterns has often attracted scientists, a mathematical or physical explanation for their common occurrence in nature is yet to be discovered ( [3]).

Spiral like gravity waves can also exist in the atmosphere. For example, as has been discussed in [4], recent observational studies confirm the existence of spiral gravity waves radiating horizontally outward from tropical cyclones. These waves are wrapped into spirals by the tangential wind of the cyclone and are described as spiral gravity waves.

In this paper, we would like to investigate if there is some correlation between atmospheric equatorial waves and the Fibonacci spiral. Equatorial waves are a key part of the tropical climate system and play a vital role in the planet’s weather and climate by transferring heat towards the poles, and cold air towards the equator. Usually such waves are associated with large-scale perturbations of the atmospheric motion extending coherently around a full longitude circle [5] [6] [7] [8] [9] . For example, Jupiter’s North Equatorial Belt (NEB) is one of the broadest and darkest belts on the planet. It is almost always a scene of notable weather formations and activity, as shown in Figure 2. Bounded by the retrograde jet at 17˚N on its northern edge, and the very fast prograde jet at 7˚N on its southern edge, the visibly dark belt does not always respect these limits. There

Figure 2. Hi-res image of the planet showing typical comparable large-scale wave structure in the southern equatorial zone called the South Equatorial Disturbance (SED) and the North Equatorial Belt (NEB), as well as small superfast spots on the SEBn jet, and rifts in both NEB and SEB. Data Source: Hubble Space Telescope 2007.

is a continuous gradient of wind speed across the NEB, as shown in Figure 2. Additionally, planetary waves are considered as an important component of the long-term mean upwelling at the tropical tropopause and the planetary wave breaking in the extratropical stratosphere ( [9] [10]). However, the exact role of equatorial planetary waves in tropical upwelling has not been resolved so far.

It is also believed that planetary waves are responsible for widespread changes in the climate system, especially in the ozone variation (see e.g. [11]). The interest in planetary waves increased dramatically increase in the 1980s when the Antarctic ozone hole was discovered ( [6] [12] [13]). The ozone layer, and in turn stratospheric ozone, plays a crucial role in protecting our life on Earth by absorbing harmful solar ultraviolet radiation ( [14] [15]). Moreover, the effects of atmospheric waves on variations of the ozone layer are also essential to the control of the stratospheric temperature via atmospheric radiative heating, as those effects make El Nino to be related with the observed rapid decelerations of the stratospheric polar vortex ( [16]). Also, it was pointed out in [11] [14] [17] [18] [19] [20] [21], that El Nino is directly related with the effects of atmospheric waves on the polar vortex. It was discussed earlier in [22] that the observed breakdowns in the polar vortex are an instability caused by atmospheric waves. To improve seasonal weather forecasting it is important to understand and improve the predictions of the dynamics of the atmospheric waves (see also [23]).

It has been discussed in [7], one method of visualizing the atmospheric waves is to measure water vapor affected by Coriolis effects detected by satellites, and in order to investigate atmospheric perturbations in a form propagating wave, a novel, robust algorithm to extract ring-shape patterns from satellite and model data has been developed in [8].

In this research, we shall analyze one particular class of the exact invariant solution associated with the nonlinear shallow water equations representing planetary equatorial waves corresponding to the free boundary problem describing the non-stationary motion of an incompressible perfect fluid propagating around a solid circle. The solution of interest was obtained by means of Lie symmetries first reported in [24] . Acquaintance with group analysis is important for constructing and investigating nonlinear mathematical models of natural and engineering problems. Numerous physical phenomena can be investigated using Lie symmetries to unearth various group invariant solutions and conservation laws that provide significant physical insight into the problem. Moreover, models of natural phenomena can often be described directly in group theoretic terms. Differential equations, conservation laws, solutions to boundary value problems, and so forth can be derived from the group invariance principle.

We consider the two dimensional motion of an incompressible perfect fluid which has a free boundary $\eta$ and a solid bottom represented by a circle of radius R. The fluid is circulating around a solid circle and confined by a free boundary. For modeling purposes, we can assume that the motion of the fluid is irrotational and the pressure on a free boundary is constant. It is postulated that that the fluid depth is small compared to the radius of the circle, as shown schematically in Figure 3.

We introduce polar coordinates $x=r\cos \theta$, $y=r\sin \theta$, where $\theta$ is a polar angle, r is the distance from the origin, and we use the following notation: R is the radius of the circle representing the radius of a planet, $h\left(t,\theta \right)=h_0+\eta \left(t,\theta \right)$, where t is time, $h_0$ is a constant undisturbed level of the equatorial atmospheric layer from the center of the planet, and $\eta \left(t\mathrm{,}\theta \right)$ is the unknown level of disturbance of a free boundary. Hence, irrotational motion of a perfect fluid is confined in the domain

${\Omega }_{\eta }=\left\{\left(r,\theta \right):0\le \theta \le 2\pi ,R\le r\le R+h\left(t,\theta \right)\right\},$

Figure 3. Schematic showing an equatorial atmospheric motion.

which is bounded by a solid circular boundary $r=R$ and a perturbed free boundary $r=R+h\left(t,\theta \right)$. For the sake of simplicity, the atmospheric flow is supposed to be irrotational and pressure on a free boundary $r=R+h\left(t,\theta \right)$ is constant. It is postulated that the unperturbed level of atmospheric “depth” $h_0$ is small compared to the radius R and the homogeneous gravity field is given by the vector $g$, which is assumed to be a constant and directed to the center of the circle.

The paper is organized as follows: First we reduce the free boundary problem to a nonlinear system of shallow water approximation. Next, the invariant solutions of the shallow water system are deduced by constructing the infinite- dimensional Lie algebra composed by the operators admitted by the system. In particular, it is noted that the system is invariant under the non-uniform scaling transformation group. As the next step, we will show that the obtained invariant solution has a Fibonacci spiral-like form and has two parameters k and $t_0$. We have defined a new model hypermarameter ${\Delta }_{k}t=t-t_0$, where t is time. The question of particular interest is: can we tune the hypermarameter in order to match the exact solution to the actual Fibonacci spiral?

Main result: It was discovered that the physically relevant part of the solution matches exactly the Fibonacci spiral.

2. Shallow Water Approximation

In what follows, it is assumed that the fluid motion is potential in the domain of the motion which allows to introduce the stream function $\psi \left(t\mathrm{,}r\mathrm{,}\theta \right)$ via

$v_r=-\frac{1}{r}\frac{\partial \psi }{\partial \theta },v_{\theta }=\frac{\partial \psi }{\partial r}.$ (1)

Notice that $\psi$ is a harmonic function in ${\Omega }_{\eta }$, since we assumed that the flow is irrotational. Next, we define the average velocity $u\left(\theta \mathrm{,}t\right)$ in terms of the stream function $\psi$ as

$u\left(\theta ,t\right)=\frac{1}{h}\stackrel{R+h}{\underset{R}{{\displaystyle \int }}}{v}_{\theta }\left(r,\theta ,t\right)\text{d}r=\frac{1}{h}\psi \left(R+h,\theta ,t\right).$ (2)

In order to reduce the number of parameters, we introduce dimensionless variables:

$\theta =\hat{\theta },r=R+h_0\hat{r},h=h_0\hat{h},t=\frac{R\hat{t}}{\sqrt{g{h}_0}},\psi =h_0\sqrt{g{h}_0}\hat{\psi },u=\sqrt{g{h}_0}\hat{u}.$ (3)

We next introduce the parameter

$\epsilon =\frac{h_0}{R}.$ (4)

The hypothesis that water is “shallow” is based on the assumption that the parameter $\epsilon$ is “small”.

In particular the relation between the boundaries of the dimensional and nondimensional variables can be summarized in the table below

Hereafter the symbol “hat” is omitted.

So, in the nondimensional variables, the free boundary is given by $r=1+\epsilon \left(1+h\right)$ and thus the functions $h\left(\theta \mathrm{,}t\right)$ and $\psi \left(r\mathrm{,}\theta \mathrm{,}t\right)$ are two unknown functions whereas $\epsilon$ is a given parameter. Although shallow water theory is usually related to the case when the water depth is small relative to the wavelengths of the waves, we find it more appropriate to choose the radius R as a natural physical scale since, in the frame of the present model, we consider waves with wavelengths of the order of the radius of a planet.

Representing the stream function $\psi$ by the series expansion: $\psi =\underset{n}{{\displaystyle \sum }}{\epsilon }^n{\psi }^{\left(n\right)}$ and application of the Lagrange’s method allows one to write the boundary problem as the following nondimensional system of nonlinear shallow water equations, which are a higher-order analogue of the Su - Gardner Equations ( [25]):

$u_t+u{u}_{\theta }+h_{\theta }+\frac{\epsilon }{2}\left(3h{u}_t-u{h}_t-u^2{h}_{\theta }+2h{h}_{\theta }\right)=0,$ (6)

$h_t+u{h}_{\theta }+h{u}_{\theta }+\epsilon h{h}_t=0,$ (7)

where the subscripts denote partial derivatives in which independent variables t and $\theta$ denote the time and the polar angle, respectively, the dependent variables are the average velocity u and the level $h>0$ of the atmosphere perturbed from $h_0$, whereas $\epsilon \ll 1$ is a small parameter.

Due to the fact that $\epsilon$ is a small parameter, the terms of order $0\left(\epsilon \right)$ can be considered as a small perturbation to the zeroth order terms (unperturbed model). Our main concern is a simplified version of the model, in which the perturbations (nonlinear terms at $\epsilon$ ) are ignored and we consider the following unperturbed system. As has been shown in [26], and [27], elimination of $u_t$ and $h_t$ from the terms of the Equations (6) and (7) and ignoring the terms with $\epsilon$ leads to the following unperturbed system:

$u_t+u{u}_{\theta }+h_{\theta }=0,$ (8)

$h_t+u{h}_{\theta }+h{u}_{\theta }=0.$ (9)

It can be checked by direct substitution that one particular exact solution of the unperturbed shallow water system (8) and (9) is

$h^0=\epsilon ,u^0=-\frac{\Gamma }{2\pi \epsilon }\ln \left(1+\epsilon \right),$ (10)

where $\Gamma =\text{const}\text{.}$ is the intensity of the vortex (source) localized at the center of the planet and is related with the rotation rate of the planet (for example, the angular velocity of the earth is $\Omega =2\pi \text{rad}/\text{day}\approx 0.73\times {10}^{-4}{\text{s}}^{-1}$ ) by the equation $\Gamma =2\pi \Omega R^2$. Since the solution (10) corresponds to the constant flow with an undisturbed circular free surface, we call it a trivial solution. To a certain extent, the above ansatz (10) can be associated with the polar vortex, which represents a very powerful whirlpool swirling steadily around the planet’s poles at all times.

Understanding and predicting break downs and overall dynamic structure of the polar vortex is important for improving seasonal forecasting (see also [23]). If the undulating path of the west to east atmospheric flow generated by the polar vortex could be predicted, then weather could be predicted too, not just for a week or two, but for an entire season [28]. However, as pointed out in [29], until now understanding of such a correlation has been based on observations and statistical modeling only rather than on the knowledge of its physical foundation. Because of the periodic seasonal breakdowns of the polar vortex and the lack of continuous data for the input of such interactions, statistical approach generally cannot provide the weather patterns accurately enough ( [30] [31] [32] [33]).

Up to the present days, our knowledge of gravity wave sources and properties in the polar region is very much limited because collecting the observations is generally difficult because of harsh natural environments ( [34]).

3. Nontrivial Similarity Solution

Detailed presentations of the theory of symmetries and invariant solutions of differential equations can be found elsewhere [24] [35] - [41]. For convenience, we summarize the basic notation from calculus of Lie group analysis in the Appendix, which represents a simplified version of the overview of basic concepts of Lie symmetry groups.

A simple inspection shows that the system (8) and (9) admits the infinite- dimensional Lie algebra composed by the operators

$\begin{array}{l}{X}_1=t\frac{\partial }{\partial \theta }+\frac{\partial }{\partial u},X_2=\theta \frac{\partial }{\partial \theta }+2h\frac{\partial }{\partial h}+u\frac{\partial }{\partial u},X_3=t\frac{\partial }{\partial t}+\theta \frac{\partial }{\partial \theta },\\ X_4=\left(2\theta -6tu\right)\frac{\partial }{\partial t}+\left(6th-3t{u}^2\right)\frac{\partial }{\partial \theta }+4hu\frac{\partial }{\partial h}+\left(u^2+4h\right)\frac{\partial }{\partial u}.\end{array}$ (11)

In particular, the operator $X_1$ is admitted due to the invariance of the system (8) and (9) under the Galilean transformation group

$\bar{\theta }=\theta +t{a}_1,\bar{u}=u+a_1,$ (12)

with an arbitrary parameter $a_1$. The operator $X_2$ indicates that the system (8) and (9) is invariant under the non-uniform scaling transformation group

$\bar{\theta }=\theta {\text{e}}^{a_2}\mathrm{,}\bar{h}=h{\text{e}}^{2a_2}\mathrm{,}\bar{u}=u{\text{e}}^{a_2}\mathrm{,}$ (13)

where $a_2$ is another arbitrary parameter and the operator $X_4$ is responsible for a specific non-scaling symmetry of the one-dimensional shallow water model (8) and (9).

The term similarity solution refers to invariant solutions based on scaling transformations (12). It can be checked by direct substitution that the shallow water model (8) and (9) is invariant under the transformation (12). One calculates the invariant $J\left(t\mathrm{,}\theta \mathrm{,}u\mathrm{,}h\right)$ of the group $X_2$ by solving the first-order linear partial differential equation

$X_{2}J=0.$ (14)

The latter Equation (14) gives three functionally independent invariants corresponding to the transformation (12):

$J_1=t,J_2=\frac{u}{\theta },J_3=\frac{h}{{\theta }^2}.$

Accordingly, we look for the invariant solution of the (8) and (9) in the form

$u=\theta U\left(t\right),h={\theta }^{2}H\left(t\right).$ (15)

Direct substitution of the presentation (15) in the system (8) and (9) yields the following nonlinear ordinary differential equations:

$\frac{\text{d}U}{\text{d}t}+U^2+2H=0,$ (16)

$\frac{\text{d}H}{\text{d}t}+3UH=0.$ (17)

In the case when $H\ne 0$, the Equation (17) can be written as the following coupled equations

$H={\text{e}}^{-3W}\mathrm{,}U=\frac{\text{d}W}{\text{d}t}\mathrm{,}$ (18)

where W satisfies the nonlinear differential equation of the second order,

$\frac{{\text{d}}^{2}W}{\text{d}{t}^2}+{\left(\frac{\text{d}W}{\text{d}t}\right)}^2+2{\text{e}}^{-3W}=0.$ (19)

Integration of the Equation (19) yields

$\frac{\text{d}W}{\text{d}t}=2{\text{e}}^{-\frac{3W}{2}}\sqrt{1+k{\text{e}}^W}\mathrm{,}$ (20)

where k is a constant and so W is given implicitly by equation

${\displaystyle \int }\frac{{\text{e}}^{3W/2}}{\sqrt{1+k{\text{e}}^W}}\text{d}W=\pm 2\left(t-t_0\right)\mathrm{,}$ (21)

where $t_0$ is an arbitrary constant.

In particular, if the function H is known, the function u can be expressed in terms of H as

$u=\pm 2\theta \sqrt{H}\sqrt{1+k{H}^{-1/3}}\mathrm{.}$ (22)

The Equation (21) provides an implicit representation of the function $W\left(t\right)$ and hence the functions $U\left(t\right)$ and $H\left(t\right)$ are (15) due to the Equations (18). We calculate the integral and distinguish the following three cases:

3.1. Case k = 0

When $k=0$ the Equation (21) has the form

$F_0\left(W\right)\mathrm{<mo>\; :\; </mo>}={\text{e}}^{\frac{3W}{2}}=\pm 3\left(t-t_0\right)\mathrm{.}$ (23)

In this case, the Equations (15) provide us with the solution

$H=\frac{1}{9{\left(t-t_0\right)}^2},U=\frac{2}{3\left(t-t_0\right)},$ (24)

and thus the exact solution of the system (8) and (9) is given by

$u=\frac{2\theta }{3\left(t-t_0\right)},h=\frac{{\theta }^2}{9{\left(t-t_0\right)}^2}\cdot$ (25)

3.2. Case k > 0

When $k>0$, the functions U and H are given by

$U=\pm 2{\text{e}}^{-\frac{3W}{2}}\sqrt{1+k{\text{e}}^W}\mathrm{,}H={\text{e}}^{-3W}\mathrm{,}$ (26)

where $k>0$ and W is given implicitly by equation

$F_{+}\left(W\right)\mathrm{<mo>\; :\; </mo>}={\text{e}}^{\frac{W}{2}}\sqrt{\frac{1}{k}+{\text{e}}^W}-\frac{1}{k}\ln \left({\text{e}}^{\frac{W}{2}}+\sqrt{\frac{1}{k}+{\text{e}}^W}\right)=\pm 2\left(t-t_0\right)\mathrm{.}$ (27)

3.3. Case k < 0

When $k<0$, the functions U and H are given by

$U=\pm 2{\text{e}}^{-\frac{3W}{2}}\sqrt{1+k{\text{e}}^W}\mathrm{,}H={\text{e}}^{-3W}\mathrm{,}$ (28)

where $k<0$ and W is given implicitly by equation

$F_{-}\left(W\right)\mathrm{<mo>\; :\; </mo>}=-{\text{e}}^{\frac{W}{2}}\sqrt{-\frac{1}{k}-{\text{e}}^W}-\frac{1}{k}\arcsin \left(\sqrt{-k}{\text{e}}^{\frac{W}{2}}\right)=\pm 2\left(t-t_0\right)\mathrm{.}$ (29)

The function $F_{-}\left(W\right)$ is defined subject to constraint

$\left|k\right|<{\text{e}}^{-W}\mathrm{.}$ (30)

For example, the function $F_{-}\left(W\right)$ is not defined for $k=0.1$.

We will approximate the functions $F_{\pm }\left(W\right)$ for large and small values of W and also for large and small values of k.

4. Asymptotic Analysis

We will analyze the asymptotic behavior of the similarity solution (15) for large and small W as well as for large and small values k.

4.1. Approximation of the Similarity Solution for Large and Small W

We next approximate the function $F_{+}\left(W\right)$ given by Equation (27) in the limit $W\to \pm \infty$. We start with $F_{+}\left(W\right)=f\left(W\right)+g\left(W\right)$, where

$f\left(W\right)={\text{e}}^{\frac{W}{2}}\sqrt{\frac{1}{k}+{\text{e}}^W}\mathrm{,}g\left(W\right)=-\frac{1}{k}\ln \left({\text{e}}^{\frac{W}{2}}+\sqrt{\frac{1}{k}+{\text{e}}^W}\right)\mathrm{.}$ (31)

First, we consider the limiting behavior of $f\left(W\right)$ :

$\underset{W\to -\infty }{\lim }f\left(W\right)=0$ (32)

and for $W\gg 0$, since $f\left(W\right)~{\text{e}}^W$, we look for $\underset{W\to +\infty }{\lim }f\left(W\right)-{\text{e}}^W$. Using $\frac{1}{u}={\text{e}}^{\frac{W}{2}}$, we find

$\underset{W\to +\infty }{\lim }f\left(W\right)-{\text{e}}^W=\underset{u\to 0}{lim}\frac{1}{u}\sqrt{\frac{1}{k}+\frac{1}{u^2}}-\frac{1}{u^2}=\underset{u\to 0}{lim}\frac{u\sqrt{\frac{u^2}{k}+1}-u}{u^3}\mathrm{.}$ (33)

After applying the L’Hopital’s rule 3 times, we get

$\underset{W\to +\infty }{\lim }f\left(W\right)=\frac{1}{2k}\mathrm{.}$ (34)

Next, with the limiting behavior of $g\left(W\right)$ :

$\underset{W\to -\infty }{\lim }g\left(W\right)=\frac{\ln k}{2k}\mathrm{.}$ (35)

For $W\gg 0$, since $g\left(W\right)~-\frac{1}{k}\ln \left(2{\text{e}}^{\frac{W}{2}}\right)$, we get

$\underset{W\to +\infty }{\lim }g\left(W\right)+\frac{\ln 2}{k}+\frac{W}{2k}=0.$ (36)

Using this, we can find two approximations for $F_{+}\left(W\right)$ :

$W\ll \mathrm{0\; <mo>\; :\; </mo>}{F}_{+}\left(W\right)\approx \frac{\ln k}{2k}\mathrm{,}$ (37)

since

$\underset{W\to -\infty }{\lim }\left(f\left(W\right)+g\left(W\right)-\left[\frac{\ln k}{2k}\right]\right)=0.$ (38)

$W\gg \mathrm{0\; <mo>\; :\; </mo>}{F}_{+}\left(W\right)\approx {\text{e}}^W+\frac{1}{k}\left(\frac{1}{2}-\ln 2-\frac{W}{2}\right)\mathrm{,}$ (39)

since

$\underset{W\to +\infty }{\lim }\left(f\left(W\right)+g\left(W\right)-\left[{\text{e}}^W+\frac{1}{k}\left(\frac{1}{2}+\ln 2-\frac{W}{2}\right)\right]\right)=0.$ (40)

This analysis shows that there are two approximations for $F_{+}\left(W\right)$, for very large W and for very small W. In particular, $\underset{W\to -\infty }{\lim }{F}_{-}\left(W\right)=0$ but $\underset{W\to +\infty }{\lim }{F}_{-}\left(W\right)$ does not exist since $F_{-}\left(W\right)$ is subject to the constraint $W<-\ln \left(-k\right)$.

4.2. Approximation of the Similarity Solution for Large and Small k

We note that the function $F_{+}\left(W\right)$ represents a better approximation of $F_0\left(W\right)$ for larger values of $\left|k\right|$ whereas the function $F_{-}\left(W\right)$ represents a better approximation of $F_0\left(W\right)$ for smaller values of $\left|k\right|$, which is demonstrated in Figure 4. In particular, the panel (a) shows the curve $F_0$ and $F_{+}$ for $k=0.1$ and $k=10$. It can be shown that $\underset{k\to \infty }{\lim }{F}_{+}=F_0$. The panel (b) compares the curve $F_0$ and $F_{-}$ for $k=0.04$ and $k=0.020$. It can be shown that $\underset{k\to 0}{\lim }{F}_{-}=F_0$, subject to the constraint (30). Finally, the panel (c) is used to compare $F_0$ with $F_{+}$ and $F_{-}$ for $k=0.02$. For this reason we call the exact solution (15) an attractor solution in the sense that $F_0$ is defined as the smallest unit which cannot be itself decomposed into two or more attractors with distinct basins of attraction. This restriction is necessary since, in general, a dynamical system may have multiple attractors, each with its own basin of attraction.

5. Nonlinear Analysis

We start by analyzing the exact attractor solution $h\left(\theta \mathrm{,}t\right)$ representing the deviation of the free boundary (see Equation (15)) form the unperturbed boundary $r=1+\epsilon$ for the case $k=0$ at different values of time t. As shown in Figure 5, the attractor solution is a decreasing function of t for a fixed value of $\theta$ but it is an increasing function of $\theta$ for a fixed value of time t. We will also apply the numerical technique to compare the attractor solution

$h=\frac{{\theta }^2}{9{\left(t-t_0\right)}^2}$ (41)

with the exact solution

$h\left(\theta \mathrm{,}t\right)={\theta }^2{\text{e}}^{-3W\left(t\right)}$ (42)

for both cases when $k>0$ and $k<0$. The results are shown in Figure 6 and Figure 7.

Figure 4. Approximation of $F_0\left(W\right)$ by $F_{+}\left(W\right)$ and $F_{-}\left(W\right)$ for larger and small values of $\left|k\right|$. We observe that $\underset{k\to \infty }{\lim }{F}_{+}=F_0$, as illustrated in panel (a) and $\underset{k\to 0}{\lim }{F}_{-}=F_0$, as illustrated in panel (b). Panel (c) is used to compare the curves $F_{+}$ and $F_{-}$ with $F_0$ for the value of $k=0.02$.

Figure 5. Attractor solution $h=\frac{{\theta }^2}{9{\left(t-t_0\right)}^2}$ for the case $k=0$ and different values of time $t=1,t=2,t=5$ and $t=10$.

Figure 6. Comparison of the solution $h\left(\theta \mathrm{,}t\right)$ for $k>0$ with the attractor solution $h=\frac{{\theta }^2}{9{\left(t-t_0\right)}^2}$ at $t=1$ and different values of k. Panel (a) shows the attractor solution (dotted red line) and the solution $h\left(\theta \mathrm{,}t\right)$ evaluated at $k=10$ (dashed blue line) and $h\left(\theta \mathrm{,}t\right)$ evaluated at $k=1$ (solid green line) and panel (a) compares these solutions for the $k=10$ and 1 as well, but evaluated at $t=2$.

Figure 7. Comparison of the solution $h\left(\theta \mathrm{,}t\right)$ for $k<0$ with the attractor solution $h=\frac{{\theta }^2}{9{\left(t-t_0\right)}^2}$ at $t=1$ and different values of k. The attractor solution is plotted by the red line and the solution $h\left(\theta \mathrm{,}t\right)$ evaluated at $k=-0.4$ is shown by the blue line whereas $h\left(\theta \mathrm{,}t\right)$ evaluated at $k=-0.1$ is shown by the green line.

In particular, Figure 6 is used to compare the solution $h\left(\theta \mathrm{,}t\right)$ given by (42) for $k>0$ with the attractor solution given by (41) at $t=1$ and different values of k. Panel (a) shows the attractor solution (dotted red line) and the solution $h\left(\theta \mathrm{,}t\right)$ given by (41) evaluated at $k=10$ (dashed blue line) and $k=1$ (solid green line). We can observe that the attractor solution represents a better approximation of the exact solution given by (41) at $k=10$ then at $k=1$, as expected from the asymptotic analysis showing that $\underset{k\to \infty }{\lim }{F}_{+}=F_0$. A similar conclusion holds for larger values of time t, as illustrated in Panel (b) of this figure, showing the solutions evaluated at $t=2$.

Similarly, Figure 7 is used to compare the nonlinear solution $h\left(\theta \mathrm{,}t\right)={\theta }^2{\text{e}}^{-3W\left(t\mathrm{,}k\right)}$ given by (42) for $k<0$ with the attractor solution

$h=\frac{{\theta }^2}{9{\left(t-t_0\right)}^2}$ given by (41) at $t=1$ and different values of k. In this plot, the

attractor solution is plotted by the red line and the solution $h\left(\theta \mathrm{,}t\right)$ evaluated at $k=-0.4$ is shown by the blue line whereas $h\left(\theta \mathrm{,}t\right)$ evaluated at $k=-0.1$ is shown by the green line. We note that the similarity solution $h\left(\theta \mathrm{,}t\right)$ is defined

for $k<0$ and $0<t<-\frac{\pi }{4k}$, which covers the whole range of $F_{-}\left(W\right)$.

6. Fibonacci Spirals

We associate the parameter k with the Fibonacci sequence $\left\{k_n\right\}$ defined by $k_n=k_{n-1}-k_{n-2}$ with $k_0=1$ and $k_1=1$. Here each value $k_n$ corresponds to

the circular arc $\theta \in \left[\frac{m\pi }{2}\mathrm{,}\frac{\left(m+1\right)\pi }{2}\right]$, where $m=0,1,2,\cdots$ so that $k_1=1$ corresponds to the circular arc $\theta \in \left[\mathrm{0,}\frac{\pi }{2}\right]$ with $m=0$. Also $k_2=1$ corresponds to the circular arc $\theta \in \left[\frac{\pi }{2}\mathrm{,}\pi \right]$ with $m=1$ and $k_3=2$ corresponds to the circular arc $\theta \in \left[\pi \mathrm{,}\frac{3\pi }{2}\right]$ with $m=2$ and so on. Next, the parameter $t_0$ is chosen for each value of k such that the endpoints of the arc segments

$\theta \in \left[\frac{m\pi }{2}\mathrm{,}\frac{\left(m+1\right)\pi }{2}\right]$ connect continuously with the next circular segment $\theta \in \left[\frac{\left(m+1\right)\pi }{2}\mathrm{,}\frac{\left(m+2\right)\pi }{2}\right]$, for any given value of time t, so that $t-t_0={\Delta }_{k}t$,

where we call ${\Delta }_{k}t$ model hypermarameter and it is tuned for the given predictive model represented by the Fibonacci spiral. We are tuning the hypermarameter in order to discover how the parameters $k_n$ and $t_0$ of the model result in the given prediction.

For example, for $t=1$, the values of ${\Delta }_{k}t$ corresponding to the first seven terms of Fibonacci sequence $\left\{k_n\right\}$ are approximated by the values shown in Table (43).

(43)

The values of ${\Delta }_{k}t$ are also shown in Figure 8 and the corresponding few circular arcs are shown in Figure 9. We observe that

${\underset{k\to \infty }{\lim }{\Delta }_k|}_{t=1}=-0.41,$ (44)

which is also seen from the Figure 8.

In particular, based on the results for the hyperparameter ${\Delta }_{k}t$ as shown in Table (43), we visualize our exact solution (15) for $h\left(\theta \mathrm{,}t\right)={\theta }^2{\text{e}}^{-3W\left(t\right)}$ for $t=1$ in Figure 9 as the sequence of circular arcs plotted at the values of time given by $t+{\Delta }_{k}t$. For example, when $k=2$, we evaluate our exact solution as a function of $\theta$ at the value of time $1+0.22=1.22$. Similarly, when $k=3$, we evaluate our exact solution as a function of $\theta$ at the value of time $1+0.287=1.287$ and so forth.

Figure 8. The values of ${\Delta }_{k}t$ calculated for the first fourteen terms of the Fibonacci sequence.

Figure 9. First few circular segments connected with each other for the first few terms of the Fibonacci sequence. The smoothness of connected is designed by choosing the corresponding values of the hypermarameter ${\Delta }_{k}t$.

Figure 10 is used to compare the attractor solution $h=\frac{{\theta }^2}{9{\left(t-t_0\right)}^2}$ with the numerical approximation of the perturbation $h\left(\theta \mathrm{,}t\right)={\theta }^2{\text{e}}^{-3W\left(t\right)}$ for $k=1$ and $\theta =\frac{\pi }{8}$ and $\theta =\frac{\pi }{4}$ as the function of time t. We note that $h\left(\theta \mathrm{,}t\right)$ represents

the deviation of the free boundary from the unperturbed state $r=1+\epsilon$, so the feasible domain of h is given by

Figure 10. Comparison of the attractor solution with $k=0$ with the numerical approximation of the perturbation $h\left(\theta \mathrm{,}t\right)$ for $k=1$ and $\theta =\frac{\pi }{8}$ and $\theta =\frac{\pi }{4}$.

${\Phi }_h=\left\{h:0\le h\le \epsilon \right\},$ (45)

where $\epsilon <1$ is a small parameter since it is assumed that the unperturbed level of atmospheric “depth” is small compared to the radius of the Earth. As Figure 10 shows, h is a decreasing function of time t, but it is an increasing function of a polar angle $\theta$, which means that our nontrivial solution is valid for only very small values of $\theta$.

Figure 11 is used to compare the free boundary given by $\eta =1+h\left(\theta ,t\right)$ and the Fibonacci spiral for the first four values of k, namely, for $k=1,1,2,3$. According to the values chosen as shown in Table (43), panel (a) shows the segment of the Fibonacci spiral (green circles) and the free boundary (blue solid line) evaluated at the values of $k=1$, $t=1$ and $t_0=0$. Panel (b) compares the segment of the Fibonacci spiral and the free boundary evaluated at the values of $k=1$, $t=1$ and $t_0=0$. Panel (c) compares the segment of the Fibonacci spiral and the free boundary evaluated at the values of $k=2$, $t=1.22$ and $t_0=-0.22$. Finally, panel (d) compares the segment of the Fibonacci spiral and the free boundary evaluated at the values of $k=3$, $t=1.287$ and $t_0=-0.287$.

As it was discussed earlier, the solution for the free boundary perturbation $h\left(\theta \mathrm{,}t\right)$ is only valid in the domain ${\Phi }_h$ defined by (45), i.e. in a very narrow band $\left|h\right|\le \epsilon <1$. For example, if $\epsilon =0.15$, we can see from Figure 10 that our

solution is valid for very small values of $\theta$ bounded by $\theta \in \left(\mathrm{0,}\frac{\pi }{4}\right)$. Thus, for example, if $\theta =\frac{\pi }{4}$, the solution is bounded by $\left|h\right|\le 0.15$. So, in this domain,

the Fibonacci spiral corresponding to $k=1$ is a very good approximation of the nontrivial solution $\eta$, as shown in panel (a) of Figure 11. For larger values of k, we observe that the solution $\eta$ diverges from the Fibonacci spiral.

Figure 11. Comparison of the free boundary $\eta =1+\epsilon h\left(\theta ,t\right)$ and the Fibonacci spiral for the first three values of k.

7. Discussions

In this paper, we have analyzed and visualized the exact invariant solution of the nonlinear simplified version of the shallow water Equations (8) and (9)

$\begin{array}{l}{u}_t+u{u}_{\theta }+h_{\theta }=0,\\ h_t+u{h}_{\theta }+h{u}_{\theta }=0,\end{array}$ (46)

which are being used to simulate equatorial atmospheric waves of planetary scales. Our model is represented by the Cauchy-Poisson free boundary problem on the nonstationary motion of a perfect incompressible fluid circulating around a vortex field approximated by a circle of a large radius and the gravity is directed to the center of the circle. The solution for the free boundary perturbation $h\left(\theta \mathrm{,}t\right)$ is only valid in a very narrow band $\left|h\right|\le \epsilon <1$, as also shown schematically in Figure 12. We have shown that within this band, the Fibonacci spiral corresponding to $k=1$ is a very good approximation of the nontrivial solution $\eta$, as shown in Figure 13. However, we observe that the solution $\eta$ diverges from the Fibonacci spiral for increasing values of k. In other words, the physically relevant part of the solution matches exactly the Fibonacci spiral.

We also remark that the higher-order shallow water equations can be derived for our free boundary model and the higher-order approximation has the following form:

$u_t+u{u}_{\theta }+h_{\theta }-\epsilon h\left(u{u}_{\theta }+\frac{h_{\theta }}{2}\right)=0,$ (47)

Figure 12. The model geometry.

Figure 13. Comparison of the free boundary given by $h=1+\epsilon \eta \left(\theta ,t\right)$ and the Fibonacci spiral evaluated at the values of $k=1$, $t=1$ and $t_0=0$.

$h_t+u{h}_{\theta }+h{u}_{\theta }-\epsilon h\left(u{h}_{\theta }+h{u}_{\theta }\right)=0.$ (48)

It is one of our goals for further studies to study the higher-order approximations (47) and (48) and see if the Fibonacci spiral can be a better approximation of the free boundary model for larger values of k.

We also note that the model (8) and (9) is invariant under the group of transformations

$X^{\infty }={\xi }^1\left(h\mathrm{,}u\right)\frac{\partial }{\partial t}+{\xi }^2\left(h\mathrm{,}u\right)\frac{\partial }{\partial \theta }\mathrm{,}$ (49)

where the functions ${\xi }^1\left(h\mathrm{,}u\right)\mathrm{,}{\xi }^2\left(h\mathrm{,}u\right)$ solve the system of first-order equations

${\xi }_h^2+{\xi }_u^1-u{\xi }_h^1=0,{\xi }_u^2-u{\xi }_u^1+h{\xi }_h^1=0.$ (50)

If we write the operator (49) as an infinitesimal transformation

$\bar{t}\approx t+{\xi }^1\left(h\mathrm{,}u\right)\mathrm{,}\bar{\theta }\approx \theta +{\xi }^2\left(h\mathrm{,}u\right)\mathrm{.}$ (51)

The Equations (50) show that the infinitesimal transformation (51) changes the variables $\left(t\mathrm{,}\theta \right)$ by adding to them an arbitrary solution of the system of linear differential equations

${\theta }_h+t_u-u{t}_h=0,{\theta }_u+h{t}_h-u{t}_u=0,$ (52)

and that the system (52) admits the infinitesimal transformation (51).

Hence, the operator (49) is admitted both by the nonlinear system (8) and (9) and by the linear system (52). This fact predicts a possibility to map the nonlinear Equations (8) and (9) to the linear system (52) by an appropriate change of variables. It is important to note that the linear system is homogeneous, i.e. invariant under the uniform dilation

$\bar{t}=t{\text{e}}^{a_5}\mathrm{,}\bar{\theta }=\theta {\text{e}}^{a_5}$

produced by the operator $X_3$ from (11).

We conclude that the operator (49) and the operator $X_3$ from (11) are responsible for the possibility of mapping of the nonlinear system (8) and (9) to the linear homogeneous system (52). Mapping the nonlinear model (8) and (9) to a linear system will be studies in the forthcoming paper.

Acknowledgements

We would like to thank the Dean of STEM & Social Sciences at Wenatchee Valley College, Holly Bringman and Electronic Media Producer Lead, Brett McGinnis for providing us with all the necessary facilities for the research, great support, and encouragement.

Appendix: Outline of Methods from Lie Group Analysis

Basic concepts from Lie group analysis of differential equations that are used in the present paper are assembled here. For further information regarding Lie groups and their applications to the theory of differential equations, the reader should consult the various classical and modern texts in the field, such as [24] [35] - [41] .

Definition of one-parameter groups. Let

${\bar{z}}^i=f^i\left(z,a\right),i=1,\cdots ,N,$ (A1)

be a one-parameter family of invertible transformations of points $z=\left(z^1,\cdots ,z^N\right)\in {ℝ}^N$ into points $\bar{z}=\left({\bar{z}}^1\mathrm{,}\cdots \mathrm{,}{\bar{z}}^N\right)\in {ℝ}^N$. Here a is a real parameter from a neighborhood of $a=0$, and we impose the condition that Transformation (A1) is an identity if and only if $a=0$, i.e.,

$f^i\left(z,0\right)=z^i,i=1,\cdots ,N.$ (A2)

The set G of transformations (A1) satisfying Condition (A2) is called a (local) one-parameter group of transformations in ${ℝ}^N$ if the successive action of two transformations is identical to the action of a third transformation from G, i.e., if the function $f=\left(f^1\mathrm{,}\cdots \mathrm{,}{f}^N\right)$ satisfies the following group property:

$f^i\left(f\left(z\mathrm{,}a\right)\mathrm{,}b\right)=f^i\left(z\mathrm{,}c\right)\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}N\mathrm{,}$ (A3)

where

$c=\phi \left(a,b\right)$ (A4)

with a smooth function $\phi \left(a\mathrm{,}b\right)$ defined for sufficiently small a and b. The group parameter a in the transformation (A1) can be changed so that the function (A4) becomes $c=a+b$. In other words, the group property (A3) can be written, upon choosing an appropriate parameter a (called a canonical parameter) in the form

$f^i\left(f\left(z\mathrm{,}a\right)\mathrm{,}b\right)=f^i\left(z\mathrm{,}a+b\right)\mathrm{.}$ (A5)

Group Generator. Let G be a group of transformations (A1) satisfying the condition (A2) and the group property (A5). Expanding the functions $f^i\left(z\mathrm{,}a\right)$ into Taylor series near $a=0$ and keeping only the linear terms in a, one obtains the infinitesimal transformation of the group G:

${\bar{z}}^i\approx z^i+a{\xi }^i\left(z\right)\mathrm{,}$ (A6)

where

${\xi }^i\left(z\right)={\frac{\partial f^i\left(z,a\right)}{\partial a}|}_{a=0},i=1,\cdots ,N.$ (A7)

The first-order linear differential operator

$X={\xi }^i\left(z\right)\frac{\partial }{\partial z^i}$ (A8)

is known as the generator of the group G.

Invariants. A function $J\left(z\right)$ is said to be an invariant of the group G if for each point $z=\left(z^1\mathrm{,}\cdots \mathrm{,}{z}^N\right)\in {ℝ}^N$ is constant along the trajectory determined by the totality of transformed points $\bar{z}$ : $J\left(\bar{z}\right)=J\left(z\right)$.

The function $J\left(z\right)$ is an invariant of the group G with Generator (A8) if and only if

$X\left(J\right)\equiv {\xi }^i\left(z\right)\frac{\partial J}{\partial z^i}=0.$ (A9)

Hence any one-parameter group has exactly $N-1$ functionally independent invariants (basis of invariants). One can take them to be the left-hand sides of $N-1$ first integrals $J_1\left(z\right)=C_1,\cdots ,J_{N-1}\left(z\right)=C_{N-1}$ of the characteristic equations for linear partial differential Equation (A9). Then any other invariant is a function of $J_1\left(z\right)\mathrm{,}\cdots \mathrm{,}{J}_{N-1}\left(z\right)$.

Invariant equations. We say that a system of equations

$F_k\left(z\right)=0,k=1,\cdots ,s$ (A10)

is invariant with respect to the group G (or admits the group G) if the transformations (A1) of the group G map any solution of Equations (A10) into a solution of the same equations, i.e.,

$F_k\left(\bar{z}\right)=0,k=1,\cdots ,s$ (A11)

whenever z solves Equations (A10). The group G with the generator (A8) is admitted by Equations (A10) if and only if

${X\left(F_k\right)|}_{\left(\text{A}10\right)}=0,k=1,\cdots ,s,$ (A12)

where the symbol |_(A10) means evaluated on the solutions of Equations (A10).

If z is a collection of independent variables $x=\left(x^1,\cdots ,x^n\right)$, dependent variables $u=\left(u^1,\cdots ,u^m\right)$ and partial derivatives $u_{\left(1\right)}=\left\{u_i^{\alpha }\right\},u_{\left(2\right)}=\left\{u_{ij}^{\alpha }\right\},\cdots$, of u with respect to x up to certain order, where

$u_i^{\alpha }=\frac{\partial u^{\alpha }}{\partial x^i},u_{ij}^{\alpha }=\frac{{\partial }^2{u}^{\alpha }}{\partial x^i\partial x^j},\cdots$

then (A10) is a system of partial differential equations

$F_k\left(x,u,u_{\left(1\right)},\cdots \right)=0,k=1,\cdots ,s.$ (A13)

Furthermore, if the transformations (A1) are obtained by the transformations of the independent and dependent variables

$\bar{x}=f\left(x,u,a\right),\bar{u}=g\left(x,u,a\right)$ (A14)

and the extension of (A14) to all derivatives $u_{\left(1\right)}$, etc. involved in the differential Equations (A13), then Equations (A11) define a group G of transformations (A14) admitted by the differential Equations (A13). In other words, an admitted group does not change the form of the system of differential Equations (A13). The generator of the admitted group G is termed an infinitesimal symmetry (or simply symmetry) of the differential Equations (A13). Equations (A12) serve for obtaining the infinitesimal symmetries and are known as the determining equations. These equations are linear and homogeneous and therefore the set L of its solutions is a vector space. Integration of determining equations often provides several linearly independent infinitesimal symmetries. Moreover, the determining equations have a specific property that guarantees that the set L is closed with respect to the commutator $\left[X_1,X_2\right]=X_1{X}_2-X_2{X}_1$. Due to this property L is called a Lie algebra. If the dimension of the vector space L is equal to r, the space is denoted by $L_r$ and is called an r-dimensional Lie algebra. An r- dimensional Lie algebra $L_r$ generates a group depending on r parameters which is called an r-parameter group.

Invariant solutions. Let the differential Equations (A13) admit a multi- parameter group G, and let H be a subgroup of G. A solution

$u^{\alpha }=h^{\alpha }\left(x\right),\alpha =1,\cdots ,m$ (A15)

of Equations (A13) is called an H-invariant solution (termed for brevity an invariant solution) if Equations (A15) are invariant with respect to the subgroup H. If H is a one-parameter group and has the generator X, then the H-invariant solutions are constructed by calculating a basis of invariants $J_1\mathrm{,}{J}_2\mathrm{,}\cdots$.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Winfree, A. (1972) Spiral Waves of Chemical Activity. Science, 175, 634-636. https://doi.org/10.1126/science.175.4022.634
[2]	Vanag, V. and Epstein, I. (2003) Segmented Spiral Waves in a Reaction-Diffusion System. PNAS, 100, 14635-14638. https://doi.org/10.1073/pnas.2534816100
[3]	Li, C., Ji, A. and Cao, Z. (2007) Stressed Fibonacci Spiral Patterns of Definite Chirality. Applied Physics Letters, 90, Article ID: 164102. https://doi.org/10.1063/1.2728578
[4]	Nolan, D. and Matthew, J. (2022) The Representation of Spiral Gravity Waves in a Mesoscale Model With Increasing Horizontal and Vertical Resolution. Journal of Advances in Modeling Earth Systems, 14, e2022MS002989. https://doi.org/10.1029/2022MS002989
[5]	Ibragimov, R., Mohazzabi, P., Roembke, R. and Van Ee, J. (2018) Asymptotic Stability of the Polar Vortex Perturbed by Harmonic Waves Describing Atmospheric Gravity Waves Circulating in an Equatorial Plane of a Spherical Planet. Mathematical Modelling of Natural Phenomena, 13, 36. https://doi.org/10.1051/mmnp/2018017
[6]	Newman, P. (2010) Chemistry and Dynamics of the Antarctic Ozone Hole. Geophysical Monograph Series, 190, 157-171. https://doi.org/10.1029/2009GM000873
[7]	Biebricher, A. (2015) Planetary Waves: The Biggest Coherent Pressure Structures in the Earth’s Atmosphere. Space for Science.
[8]	Gong, J., Yue, J. and Wu, D. (2015) Global Survey of Concentric Gravity Waves in AIRS Images and ECMWF Analysis. GJGR Atmospheres, 120, 2210-2228. https://doi.org/10.1002/2014JD022527
[9]	Holton, J., Haynes, P., McIntyre, M., Douglass, A., Rood, C. and Pfister, L. (1995) Stratosphere-Troposphere Exchange. Reviews of Geophysics, 33, 403-439. https://doi.org/10.1029/95RG02097
[10]	Yulaeva, E. and Wallace, J. (1994) The Signature of ENSO in Global Temperatures and Precipitation Fields Derived from the Microwave Sounding Unit. Journal of Climate, 7, 1719-1736. https://doi.org/10.1175/1520-0442(1994)007<1719:TSOEIG>2.0.CO;2
[11]	Zhang, Y., Li, J. and Zhou, J. (2017) The Relationship between Polar Vortex and Ozone Depletion in the Antarctic Stratosphere during the Period 1979-2016. Advances in Meteorology, 2017, Article ID: 3078079. https://doi.org/10.1155/2017/3078079
[12]	Newman, P. and Schoeberl, M. (2003) Middle Atmosphere: Polar Vortex. In: Holton, J.R., Pyle, J. and Curry, J.A., Eds., Encyclopedia of Atmospheric Sciences, Academic, San Diego, 1321-1328. https://doi.org/10.1016/B0-12-227090-8/00228-1
[13]	Shindell, D.T. and Schmidt, G.A. (2004) Southern Hemisphere Climate Response to Ozone Changes and Greenhouse Gas Increases. Geophysical Research Letters, 31, L18209. https://doi.org/10.1029/2004GL020724
[14]	Lubin, D. and Jensen, E. (2002) Effects of Clouds and Stratospheric Ozone Depletion on Ultraviolet Radiation Trends. Nature, 377, 710-713. https://doi.org/10.1038/377710a0
[15]	Chipperfield, M., Dhomse, S., Feng, W., McKenzie, R., Velders, G. and Pyle, J. (2015) Quantifying the Ozone and Ultraviolet Benefits Already Achieved by the Montreal Protocol. Nature Communications, 6, 7233. https://doi.org/10.1038/ncomms8233
[16]	Baldwin, M. and Dunkerton, T. (2001) Stratospheric Harbingers of Anomalous Weather Regimes. Science, 80, 581-584. https://doi.org/10.1126/science.1063315
[17]	El Niño and La Niña (2016) New Zealand’s National Institute of Water and Atmospheric Research.
[18]	Becker, J. (2016) How Much Do El Niño and La Niña Affect Our Weather? This Fickle and Influential Climate Pattern Often Gets Blamed for Extreme Weather. A Closer Look at the Most Recent Cycle Shows That the Truth Is More Subtle. Scientific American.
[19]	Angell, J. and Korshover, J. (1977) Variation in Size and Location of the 300 mb North Circumpolar Vortex between 1963 and 1975. Monthly Weather Review, 105, 19-25. https://doi.org/10.1175/1520-0493(1977)105<0019:VISALO>2.0.CO;2
[20]	Sun, L., Wu, H. and Li, X. (2006) Our Understanding of Arctic Vortex. Chinese Journal of Polar Research, 18, 52-62.
[21]	McCreary, J. (1976) Eastern Tropical Ocean Response to Changing Wind Systems: With Applications to El Niño. Journal of Physical Oceanography, 6, 632-545. https://doi.org/10.1175/1520-0485(1976)006<0632:ETORTC>2.0.CO;2
[22]	Murray, F. (1960) Dynamic Stability in the Stratosphere. Journal of Geophysical Research, 65, 3273-3305. https://doi.org/10.1029/JZ065i010p03273
[23]	Alexander, M., Tsuda, T. and Vincent, R. (2002) Latitudinal Variations Observed in Gravity Waves with Short Vertical Wavelengths. Journal of the Atmospheric Sciences, 59, 1394-1404. https://doi.org/10.1175/1520-0469(2002)059<1394:LVOIGW>2.0.CO;2
[24]	Lie, S. (1891) Vorlesungen uber Differentialgleichungen Mit Bekannten Infinitesimalen Transformationen, Bearbeited und herausgegeben von Dr. G. Scheffers, B.G. Teubner, Leizig.
[25]	Su, C. and Gardner, C. (1968) Derivation of the Korteweg-de Vries Equations and Burgers Equations. Journal of Mathematical Physics, 10, 536-539. https://doi.org/10.1063/1.1664873
[26]	Ibragimov, R. (2011) Nonlinear Viscous Fluid Patterns in a Thin Rotating Spherical Domain and Applications. Physics of Fluids, 23, Article ID: 123102. https://doi.org/10.1063/1.3665132
[27]	Ibragimov, R. and Gunag, L. (2015) Splitting Phenomenon of a Higher-Order Shallow Water Theory Associated with a Longitudinal Planetary Waves. Dynamics of Atmospheres and Oceans, 69, 1-11. https://doi.org/10.1016/j.dynatmoce.2014.10.003
[28]	Stanford University (2018) Pulse of the Polar Vortex Revealed: A Key to Mapping Future Storms. ScienceDaily.
[29]	Sheshadri, A., Plumb, R., Lindgren, E. and Domeisen, D. (2018) The Vertical Structure of Annular Modes. Journal of the Atmospheric Sciences, 75, 3507-3519. https://doi.org/10.1175/JAS-D-17-0399.1
[30]	Yang, X. (2016) Variability of the Northern Circumpolar Vortex and Its Association with Climate Anomaly in China. Acta Meteorologica Sinica, 26, 135-142.
[31]	Sigmond, M., Scinocca, J., Kharin, V. and Shepherd, T. (2013) Enhanced Seasonal Forecast Skill Following Stratospheric Sudden Warmings. Nature Geoscience, 6, 98-102. https://doi.org/10.1038/ngeo1698
[32]	Shi, C., Xu, T., Guo, D. and Pan, Z. (2017) Modulating Effects of Planetary Wave 3 on a Stratospheric Sudden Warming Event in 2005. Journal of the Atmospheric Sciences, 74, 1549-1559. https://doi.org/10.1175/JAS-D-16-0065.1
[33]	Shi, C., Xu, T., Li, H. and Gao, Y. (2017) The Role of Rossby-Wave Propagation in a North American Extreme Cold Event. Advances in Meteorology, 2017, Article ID: 4635849. https://doi.org/10.1155/2017/4635849
[34]	Sato, K. and Yoshiki, M. (2008) Gravity Wave Generation around the Polar Vortex in the Stratosphere Revealed by 3-Hourly Radiosonde Observations at Syowa Station. Journal of the Atmospheric Sciences, 65, 3719. https://doi.org/10.1175/2008JAS2539.1
[35]	Cohen, A. (1911) An Introduction to the Lie Theory of One-Parameter Groups with Applications to the Solution of Differential Equations. D.C. Heath, New York.
[36]	Dickson, L.E. (1924) Differential Equations from the Group Standpoint. Annals of Mathematics, 25, 287. https://doi.org/10.2307/1967773
[37]	Ovsyannikov, L.V. (1962) Group Properties of Differential Equations. Siberian Branch, USSR Academy of Sciences, Novosibirsk.
[38]	Olver, P.J. (1986) Applications of Lie Groups to Differential Equations. Springer- Verlag, New York. https://doi.org/10.1007/978-1-4684-0274-2
[39]	Bluman, G. and Kumei, S. (1987) On Invariance Properties of the Wave Equations. Journal of Mathematical Physics, 28, 307. https://doi.org/10.1063/1.527659
[40]	Stephani, H. (1989) Differential Equations: Their Solution Using Symmetries. Cambridge Univ. Press, Cambridge. https://doi.org/10.1017/CBO9780511599941
[41]	Kovalev, V., Pustovalov, V. and Senashov, S. (1993) Lie-Bäcklund Group for the Nonlinear Geometrical Optics Equations. Differential Equations, 29, 1751-1764.