Visualization of Exact Invariant Solutions Associated with Atmospheric Waves in a Thin Circular Layer

Ranis N. Ibragimov^{1}^{*},
Lauren D. Mongrain^{1},
Benjamin Stimmel^{1},
Olga Trozkaya^{2},
Guang Lin^{3},
Sheng Zhang^{3},
Vesselin Vatchev^{4},
Daniel Stankiewicz^{5}

Show more

1. Introduction

Planetary waves naturally occur in rotating fluids largely due to the rotation of the planet. These waves affect the planet’s weather and climate by transferring heat towards the poles, and cold air towards the equator. Equatorial planetary waves are a key part of the tropical climate system and they are associated with large-scale perturbations of the atmospheric motion extending coherently around a full longitude circle [1]. Equatorial 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 ( [2] [3] ). 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. [4] ). The interest in planetary waves increased dramatically in the 1980s when the Antarctic ozone hole was discovered ( [5] [6] [7] ). The ozone layer, and in turn stratospheric ozone play a crucial role in protecting our life on Earth by absorbing harmful solar ultraviolet radiation ( [8] [9] ). 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 ( [10] ). Also, it was pointed out in [11] [12] [13] [14] [8] [4] and [15], that El Nino is directly related with the effects of atmospheric waves on the polar vortex. It was discussed earlier in [16] 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 [17] ).

It has been discussed in [18],* *one method of visualizing the atmospheric waves is to measure water vapor affected by Coriolis effects detected by satellites, as shown in Figure 1. 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 [19].

In this research, we aim to visualize the set of exact invariant solutions of nonlinear shallow water equations associated with planetary equatorial waves

Figure 1. Water vapour measurements by satellite uncover a planetary wave structure. Credit: NASA. Also available at: https://spaceforscience.wordpress.com/2015/04/22/planetary-waves-the-biggest-coherent-pressure-structures-in-the-earths-atmosphere.

corresponding to the Cauchy-Poisson free boundary problem describing the nonstationary motion of an incompressible perfect fluid propagating around a solid circle.

We consider a 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 2.

We introduce polar coordinates
$x=r\mathrm{cos}\left(\theta \right)$,
$y=r\mathrm{sin}\left(\theta \right)$, 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\mathrm{,}\theta \right)={h}_{0}+\eta \left(t\mathrm{,}\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\mathrm{,}\theta \right)\mathrm{:0}\le \theta \le 2\pi \mathrm{,}R\le r\le R+h\left(t\mathrm{,}\theta \right)\right\}\mathrm{,}$

which is bounded by a solid circular boundary
$r=R$ and a perturbed free boundary
$r=R+h\left(t\mathrm{,}\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\mathrm{,}\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.

Figure 2. Schematic showing an equatorial atmospheric motion.

2. Free Boundary Model

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},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{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.

As has been shown in [20], the model describing planetary equatorial atmospheric motion around an equatorial plane can be written as the following free boundary problem:

$\frac{{\partial}^{2}\psi}{\partial {\theta}^{2}}+{r}^{2}\frac{{\partial}^{2}\psi}{\partial {r}^{2}}+r\frac{\partial \psi}{\partial r}=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(R<r<R+h\right),$ (2)

$\psi \left(R,\theta ,t\right)=0,$ (3)

$\psi \left(R+h,\theta ,t\right)=u\left(\theta ,t\right)h,$ (4)

$\frac{{\partial}^{2}\psi}{\partial t\partial r}-\frac{1}{{r}^{2}}\frac{\partial h}{\partial \theta}\frac{{\partial}^{2}\psi}{\partial t\partial \theta}+\frac{1}{2r}\frac{\partial}{\partial \theta}\left[\frac{1}{{r}^{2}}{\left(\frac{\partial \psi}{\partial \theta}\right)}^{2}+{\left(\frac{\partial \psi}{\partial r}\right)}^{2}\right]+\frac{g}{r}\frac{\partial h}{\partial \theta}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(r=R+h\right),$ (5)

$r\frac{\partial h}{\partial t}+\frac{\partial}{\partial \theta}\left(uh\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(r=R+h\right),$ (6)

where the average velocity $u\left(\theta \mathrm{,}t\right)$ is defined 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).$ (7)

Shallow Water Approximation

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

$\begin{array}{l}\theta =\stackrel{^}{\theta},\text{\hspace{0.17em}}\text{\hspace{0.17em}}r=R+{h}_{0}\stackrel{^}{r},\text{\hspace{0.17em}}\text{\hspace{0.17em}}h={h}_{0}\stackrel{^}{h},\text{\hspace{0.17em}}\text{\hspace{0.17em}}t=\frac{R\stackrel{^}{t}}{\sqrt{g{h}_{0}}},\\ \psi ={h}_{0}\sqrt{g{h}_{0}}\stackrel{^}{\psi},\text{\hspace{0.17em}}\text{\hspace{0.17em}}u=\sqrt{g{h}_{0}}\stackrel{^}{u}.\end{array}$ (8)

We next introduce the parameter

$\epsilon =\frac{{h}_{0}}{R}.$ (9)

Of course, water is shallow if the parameter
$\epsilon $ is small. So, in the present model (2)-(6), the functions
$\eta \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.

As has been shown in [21] * *that 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 reduce the Cauchy-Poisson free boundary problem (2)-(6) to the following nondimensional system of nonlinear shallow water equations, which are a higher-order analogue of the *Su - Gardner equations* ( [22] ):

${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,$ (10)

${h}_{t}+u{h}_{\theta}+h{u}_{\theta}+\epsilon h{h}_{t}=0,$ (11)

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
$\epsilon \ll 1$ of the atmosphere perturbed from
${h}_{0}$, whereas
$\epsilon \ll 1$ is a small parameter.

We eliminate ${u}_{t}$ and ${h}_{t}$ from the terms of the Equations (10)-(11) and ignore the terms with ${\epsilon}^{2}$. This leads to the following system:

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

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

Due to the fact that $\epsilon $ is a small parameter, the terms of order $0\left(\epsilon \right)$ in (12)-(13) can be considered as 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:

${u}_{t}+u{u}_{\theta}+{h}_{\theta}=0,$ (14)

${h}_{t}+u{h}_{\theta}+h{u}_{\theta}=0.$ (15)

It can be checked by direct substitution that one particular exact solution of the unperturbed shallow water system (14)-(15) is

${h}^{0}=\epsilon ,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{u}^{0}=-\frac{\Gamma}{2\pi \epsilon}\mathrm{ln}\left(1+\epsilon \right),$ (16)

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{\hspace{0.17em}}\text{rad}/\text{day}\approx 0.73\times {10}^{-4}{\text{s}}^{-1}$ ) by the equation
$\Gamma =2\pi \Omega {R}^{2}$. The solution (16) is visualized in Figure 3. In this plot, we use the following values of the parameters:
$R=637.1\times {10}^{6}\left[\text{cm}\right]$, which corresponds to the radius of the Earth,
${h}_{0}=12\times {10}^{6}$, which corresponds to thermosphere and the parameter
$\Gamma $ is normalized to one. Since the atmospheric layer is too small compared the radius of the Earth, we introduce the “*visualization parameter*
$\nu =10$ ” and do the rescaling
$\epsilon \to \nu \epsilon $.

The exact solution $\left({h}^{0}\mathrm{,}{u}^{0}\right)$ given by (16) is shown in Figure 3 and corresponds to the constant flow with an undisturbed circular free surface and the negative sign in ${u}^{0}$ corresponds to a counter-clockwise rotation. In particular, the plot shown in Figure 3 is based on the distance between an arbitrary point at

Figure 3. Visualization of the exact solution (16) with the following values of the parameters: $R=637.1\times {10}^{6}\left[\text{cm}\right]$, which corresponds to the radius of the Earth, ${h}_{0}=12\times {10}^{6}$, which corresponds to thermosphere and the parameter $\Gamma $ is chosen so that ${u}^{0}=1$ in formula (16). Since the atmospheric layer is too small compared the radius of the Earth, we introduce the “visualization parameter $\nu =10$ ” and we do the rescaling $\epsilon \to \nu \epsilon $.

the free surface *h* and the fixed origin and the angle
$\theta $, as demonstrated by Figure 2. We visualize the flow in a color scale from −1 to 10^{8}. Since
$\Gamma $ is a constant, without loss of generality we set it to make
${u}^{0}$ to be normalized to one. At certain extent, the above ansatz (16) can be associated with the polar vortex, which represents a very powerful whirlpool swirling steadily around the planet’s poles at all times. The term polar vortex is used to describe several different features in the atmosphere. It most commonly refers to a planetary-scale mid- to high-latitude circumpolar circulation (see [23] ). A polar vortex is an upper level low-pressure area lying near the Earth’s pole. There are two polar vortices in the Earth’s atmosphere, which lie over the North and South Poles. Each polar vortex is a persistent, large-scale, low pressure zone that rotates counter-clockwise at the North Pole (called a cyclone), and clockwise at the South Pole. One spectacular example of circulating atmospheric gravity waves is demonstrated in Figure 4 showing an atmospheric stream that follows a hexagon-shaped path at the north pole of Saturn (e.g. [24] ). Saturn’s hexagon is a persisting hexagonal cloud pattern around the north pole of Saturn, located at about 78˚N. The sides of the hexagon are about 13,800 km long, which is more

Figure 4. Images from Cassini, made possible only as Saturn’s north pole emerged from winter darkness, shows new details of a jet stream that follows a hexagon-shaped path and has long puzzled scientists. Image Credit: NASA/JPL-Caltech/Space Science Institute. More Cassini images can be found at http://ciclops.org.

than the diameter of Earth (about 12,700 km), which is also shown schematically in this Figure.

The hexagon was originally discovered in images taken by Voyager spacecraft in the early 1980s. Since 2006, the Cassini Visual and Infrared Mapping Spectrometer (VIMS) instrument has been observing the hexagon at infrared wavelengths, but at lower spatial resolution than these visible light images. This phenomenon showing waves that can be seen traveling along hexagon remains unexplained.

Understanding and predicting break downs and overall dynamic structure of the polar vortex is important for improving seasonal forecasting (see also [17] ). 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 [25]. However, as pointed out in [26], until now understanding of a such 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 ( [27] [28] [29] [30] ).

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 ( [31] ).

3. Exact Invariant Solutions

Computation of invariant solutions of the shallow water system (14)-(15) is based on finding the infinite-dimensional Lie algebra.^{1}

Basic concepts from Lie group analysis of differential equations used in the present paper are assembled here. If the reader wants to have more information about Lie groups and their applications in the theory of differential equations he has a big choice of classical and modern texts in this field, e.g. [32] - [41] *.*

Definition of one-parameter groups. Let

${\stackrel{\xaf}{z}}^{i}={f}^{i}\left(z,a\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=1,\cdots ,N,$ (17)

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

${f}^{i}\left(z,0\right)={z}^{i},\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=1,\cdots ,N.$ (18)

The set *G* of transformations (17) satisfying Condition (18) is called a (local) one-parameter group of transformations in
${\mathbb{R}}^{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},\cdots ,{f}^{N}\right)$ satisfies the following group property:

${f}^{i}\left(f\left(z,a\right),b\right)={f}^{i}\left(z,c\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=1,\cdots ,N,$ (19)

where

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

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

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

Group Generator. Let *G* be a group of transformations (17) satisfying the condition (18) and the group property (21). 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*:

${\stackrel{\xaf}{z}}^{i}\approx {z}^{i}+a{\xi}^{i}\left(z\right)\mathrm{,}$ (22)

where

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

The first-order linear differential operator

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

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 {\mathbb{R}}^{N}$ is constant along the trajectory determined by the totality of transformed points
$\stackrel{\xaf}{z}$ :
$J\left(\stackrel{\xaf}{z}\right)=J(\; z\; )$

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

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

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 (25). 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,\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=1,\cdots ,s$ (26)

is invariant with respect to the group *G* (or *admits* the group *G*) if the transformations (17) of the group *G* map any solution of Equation (26) into a solution of the same equations, *i.e.*,

${F}_{k}\left(\stackrel{\xaf}{z}\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=1,\cdots ,s$ (27)

whenever *z* solves Equation (26). The group *G* with the generator (24) is admitted by Equation (26) if and only if

${X\left({F}_{k}\right)|}_{\left(26\right)}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=1,\cdots ,s,$ (28)

where the symbol ${\mathrm{|}}_{\left(26\right)}$ means evaluated on the solutions of Equation (26).

If *z* is a collection of independent variables
$x=\left({x}^{1}\mathrm{,}\cdots \mathrm{,}{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}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{u}_{ij}^{\alpha}=\frac{{\partial}^{2}{u}^{\alpha}}{\partial {x}^{i}\partial {x}^{j}},\cdots $ (29)

then (26) is a system of partial differential equations

${F}_{k}\left(x,u,{u}_{\left(1\right)},\cdots \right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=1,\cdots ,s.$ (30)

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

$\stackrel{\xaf}{x}=f\left(x,u,a\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\xaf}{u}=g\left(x,u,a\right)$ (31)

and the extension of (31) to all derivatives
${u}_{\left(1\right)}$, etc. involved in the differential equations (30), then Equation (27) define a group *G* of transformations (31) admitted by the differential Equation (30). In other words, an admitted group does not change the form of the system of differential Equation (30). The generator of the admitted group *G* is termed an *infinitesimal symmetry* (or simply *symmetry*) of the differential Equation (30). Equation (28) serves for obtaining the infinitesimal symmetries and is 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}\mathrm{,}{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

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

${u}^{\alpha}={h}^{\alpha}\left(x\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\alpha =1,\cdots ,m$ (32)

of Equation (30) is called an *H*-invariant solution (termed for brevity an *invariant solution*) if Equation (32) 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 $.

3.1. Galilean Transformation

The functionally independent invariants of the operator

$X=t\frac{\partial}{\partial \theta}+\frac{\partial}{\partial u}$ (33)

provide the following particular exact solution of the system (14)-(15):

$u\left(t,\theta \right)=\frac{1}{t}\left(\theta +B\right),\text{\hspace{1em}}h\left(t,\theta \right)=\frac{A}{t}\text{\hspace{0.05em}},$ (34)

where *A* and*B* are arbitrary constants.

The exact solution (34) is visualized in Figure 5. For the sake of simplicity, the constants *A* and *B* are set to be one. Following the notation in Figure 2, the physical meaning of *u* is an averaged atmospheric flow velocity confined between the solid surface of the earth (
$r=R$ ) and atmospheric boundary
$r=h$. According to our plot shown in Figure 3, the unperturbed atmospheric layer is
$h=\epsilon $ and the exact solution for *h* in (34) is
$h=1/t$. So, in Figure 5, we choose
$t=1/\epsilon $. In this case,
$r=h$. Also, for fixed value of time *t*, function *u* depends on
$\theta $ only (no radial dependence). Correspondingly, for the fixed value of *t*, the solution (34) is visualized as
$\theta $, dependent flow confined in an unperturbed circular layer. We visualize the flow in the same color scheme as in Figure 3. The minimum value of *u* occurs at
$\theta =0$ and then *u* is monotonically increasing to its maximum value 10^{8} attained at
$\theta =2\pi $. Then, due to periodic nature of the flow, the velocity goes to the minimum value right after
$\theta =2\pi $ and increases as values of
$\theta $ increase.

3.2. Similarity Solution

The term *similarity solution* refers to invariant solutions based on scaling symmetries [42]. The system (14)-(15) admits the following infinitesimal operator

$X=\theta \frac{\partial}{\partial \theta}+2h\frac{\partial}{\partial h}+u\frac{\partial}{\partial u},$

Figure 5. Visualization of exact solution (34) with the following values of the parameters: $R=637.1\times {10}^{6}\left[\text{cm}\right]$, which corresponds to the radius of the Earth, ${h}_{0}=12\times {10}^{6}$, which corresponds to thermosphere and the parameter $\Gamma $ is chosen so that ${u}^{0}=1$ in formula (16). Again, since the atmospheric layer is too small compared the radius of the Earth, we introduce the “visualization parameter $\nu =10$ ” and we do the rescaling $\epsilon \to \nu \epsilon $.

which provides the following functionally independent invariants:

${J}_{1}=t,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{J}_{2}=\frac{u}{\theta},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{J}_{3}=\frac{h}{{\theta}^{2}}.$ (35)

Accordingly, the system (14)-(15) admits the invariant solution in the form

$u=\theta U\left(t\right),\text{\hspace{1em}}h={\theta}^{2}H\left(t\right).$ (36)

Direct substitution of the presentation (36) in the system (14)-(15) yields the following nonlinear ordinary differential equations:

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

${\left(\frac{\text{d}H}{\text{d}t}\right)}^{2}+3UH=0.$ (38)

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

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

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.$ (40)

Integration of the Equation (40) yields

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

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

$\int}\frac{{\text{e}}^{3W/2}}{\sqrt{1+k\text{\hspace{0.05em}}{\text{e}}^{W}}}\text{d}W=\pm 2\left(t-{t}_{0}\right)\mathrm{,$ (42)

where ${t}_{0}$ is an arbitrary constant.

The equation (42) provides an implicit representation of the function $W\left(t\right)$ and hence the functions $U\left(t\right)$ and $H\left(t\right)$ is (36) due to the Equation (39). We calculate the integral and distinguish the following three cases:

3.2.1. Case *k *= 0

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

${F}_{0}\left(W\right)\mathrm{:}={\text{e}}^{\frac{3W}{2}}=\pm 3\left(t-{t}_{0}\right)\mathrm{.}$ (43)

In this case, the Equation (36) provide us with the solution

$H=\frac{1}{9{\left(t-{t}_{0}\right)}^{2}}\mathrm{,}\text{\hspace{1em}}U=\frac{2}{3\left(t-{t}_{0}\right)}$ (44)

and thus the exact solution of the system (14)-(15) is given by

$u=\frac{2\theta}{3\left(t-{t}_{0}\right)}\mathrm{,}\text{\hspace{1em}}h=\frac{{\theta}^{2}}{9{\left(t-{t}_{0}\right)}^{2}}\cdot $ (45)

The exact solution (45) is visualized in Figure 6 for the fixed values of
${t}_{0}=0.1$ and
$t=1.4>{t}_{0}$. The other values of parameters and scaling are the same as in Figure 3 and Figure 5. Also, for fixed values of time *t* and *t*_{0}, functions *u* and *h* depends on
$\theta $ only. Correspondingly, for the fixed values of *t* and *t*_{0}, the solution (34) is visualized as
$\theta $ -dependent flow confined in an
$\theta $ -dependent unperturbed circular layer. We visualize the flow in the same color scheme as in Figure 3. The minimum value of the outer atmospheric level
$h\left(\theta \right)$ occurs at
$\theta =0$ and then *h* is monotonically increasing to some maximum value at
$\theta =2\pi $. Then, due to periodic nature of the flow, the atmospheric level *h* goes to the minimum value right after
$\theta =2\pi $ and increases as values of
$\theta $ increase. The behavior of the velocity *u* remains the same is shown in Figure 5.

Figure 7 is used to show the snapshots of the exact solution (45) using the same values of parameters and scaling as in Figure 6 but at different times ${t}_{0}=0.1$ (top left), ${t}_{1}=1$ (top right), ${t}_{2}=1.1$ (bottom left) and ${t}_{3}=1.3$ (bottom right).

We observe the same qualitative behavior of *u* and *h* as in Figure 6. However, we note two major differences due to time variation. First, as *t* increases, the level of atmospheric elevation above the Earth surface decreases at
$\theta =0$. Second, the value of average velocity *u* at fixed values of
$\theta $ becomes smaller at larger values of time *t.*

Figure 6. Visualization of the exact solution (45) for the fixed value $t=1.4$. The values of parameters $R\mathrm{,}{h}_{0}\mathrm{,}\Gamma $ and $\nu $ are the same as in Figure 3 and Figure 5.

Figure 7. Snapshots of the exact solution (45) using the same values of parameters and scaling as in Figure 6 but at different times ${t}_{0}=0.1$ (top left), ${t}_{1}=1$ (top right), ${t}_{2}=1.1$ (bottom left) and ${t}_{3}=1.3$ (bottom right).

3.2.2. Case *k *> 0

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

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

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

${F}_{+}\left(W\right)\mathrm{:}={\text{e}}^{\frac{W}{2}}\sqrt{\frac{1}{k}+{\text{e}}^{W}}-\frac{1}{k}\mathrm{ln}\left({\text{e}}^{\frac{W}{2}}+\sqrt{\frac{1}{k}+{\text{e}}^{W}}\right)=\pm 2\left(t-{t}_{0}\right)\mathrm{.}$ (47)

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

$u=\pm 2\theta \text{\hspace{0.05em}}\sqrt{H}\sqrt{1+k{H}^{-1/3}}\mathrm{.}$ (48)

3.2.3. Case *k *< 0

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

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

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

${F}_{-}\left(W\right)\mathrm{:}=-{\text{e}}^{\frac{W}{2}}\sqrt{-\frac{1}{k}-{\text{e}}^{W}}-\frac{1}{k}\mathrm{arcsin}\left(\sqrt{-k}{\text{e}}^{\frac{W}{2}}\right)=\pm 2\left(t-{t}_{0}\right)\mathrm{.}$ (50)

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

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

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

Approximation of the similarity solution

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 8. 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 (51). Finally, the panel (c) is used to compare ${F}_{0}$ with ${F}_{+}$ and ${F}_{-}$ for $k=0.02$.

This observation hallows us to approximate the implicit presentations (47) for
${F}_{+}$ and (50) for
${F}_{-}$ by their linear approximations in *k* in the vicinity of the

point
$k=0$ at which, according to (43),
${\text{e}}^{\frac{3{W}_{0}}{2}}=\pm 3\left(t-{t}_{0}\right)$, where
${W}_{0}={W|}_{k=0}$ is constant for some fixed values of *t*_{0} and *t*. Thus

${W}_{0}=\frac{2}{3}\mathrm{ln}{\left(t-{t}_{0}\right)}^{3},\text{\hspace{0.17em}}\text{\hspace{0.17em}}t>{t}_{0}.$ (52)

Consider the case when $k>0$. We denote

${m}_{1}=\pm \frac{2}{P},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\alpha =\frac{1}{P}\left(\mp 2{t}_{0}-{a}_{1}{\text{e}}^{\frac{{W}_{0}}{2}}\right)+{W}_{0},$ (53)

Figure 8. 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$.

where

$P=\frac{1}{4}\frac{{\text{e}}^{\frac{{W}_{0}}{2}}{a}_{1}\left(-{\text{e}}^{\frac{{W}_{0}}{2}}+\left(k{\text{e}}^{{W}_{0}}+2{k}^{2}{\text{e}}^{{W}_{0}}\right)\left[{a}_{2}+{\text{e}}^{\frac{{W}_{0}}{2}}\right]\right)}{\left(1+k{\text{e}}^{{W}_{0}}\right)\left({\text{e}}^{\frac{{W}_{0}}{2}}+{a}_{2}\right){a}_{1}\sqrt{k}}+\frac{1}{2}{\text{e}}^{\frac{{W}_{0}}{2}}{a}_{1},$ (54)

${a}_{1}=\sqrt{\frac{1+k{\text{e}}^{{W}_{0}}-\mathrm{ln}\left({\text{e}}^{\frac{{W}_{0}}{2}}+\sqrt{\frac{1+k{\text{e}}^{{W}_{0}}}{k}}\right)}{k}}$ (55)

and ${W}_{0}$ is given by (52).

Thus, expanding
${F}_{+}\left(W\right)$ in Taylor series in the vicinity of the point
$k=0$, we get the following linear approximation for *W*:

$W={m}_{1}t+\alpha ,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}k>0.$ (56)

Now we consider the case when $k<0$. We denote

${m}_{2}=\pm \frac{2}{{T}_{2}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\beta =\frac{1}{{T}_{2}}\left(\mp 2{t}_{0}+{T}_{1}\right)+{W}_{0},$ (57)

where

${T}_{1}={\text{e}}^{\frac{{W}_{0}}{2}}b+\frac{1}{k}\mathrm{arcsin}\left[\sqrt{-k}{\text{e}}^{\frac{{W}_{0}}{2}}\right],$ (58)

${T}_{2}=\frac{1}{2}\frac{{\text{e}}^{\frac{{W}_{0}}{2}}\left(1+k\left[{b}^{2}+{\text{e}}^{{W}_{0}}\right]+2{k}^{2}{b}^{2}{\text{e}}^{{W}_{0}}\right)}{{\left(-k\right)}^{\frac{1}{4}}{b}^{\frac{3}{2}}k},$ (59)

in which

$b=\sqrt{-\frac{1+k{\text{e}}^{{W}_{0}}}{k}}.$ (60)

and ${W}_{0}$ is given by (52).

Thus, expanding
${F}_{-}\left(W\right)$ in Taylor series in the vicinity of the point
$k=0$, we get the following linear approximation for *W*:

$W={m}_{2}t+\beta ,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}k<0.$ (61)

4. Conclusions

In this paper, we have visualized some exact invariant solutions of the nonlinear simplified version of the shallow water equations

${u}_{t}+u{u}_{\theta}+{h}_{\theta}=0,$ (62)

${h}_{t}+u{h}_{\theta}+h{u}_{\theta}=0,$ (63)

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. We also found some additional invariant solutions based on similarity solutions. In particular, the nontrivial forms of the motion can be represented by the following linear combination of the Galilean and scaling infinite-dimensional Lie algebra of operators:

$X=\left(2\theta -6tu\right)\frac{\partial}{\partial t}+\left(6ht-3t{u}^{2}\right)\frac{\partial}{\partial \theta}+4hu\frac{\partial}{\partial h}+\left(4h+{u}^{2}\right)\frac{\partial}{\partial u}.$ (64)

The operator (64) provides a more complicated form of *h* than that given in (36), namely for

$h=m{H}^{1/3}+{\theta}^{2}H\left(t\right)\mathrm{,}\text{\hspace{1em}}m=\text{const}\text{.}$ (65)

We will study more complicated forms of similarity solutions based on the operator (64) in the forthcoming studies.

We also remark that the exact solutions (34) and (45) represent only “local” physically relevant solutions, *i.e.* they valid for small variations of
$\theta $ and *t*. Also, these solutions are associated with the unperturbed version of the complete

Figure 9. Schematic showing multiple atmospheric waves approximating an equatorial atmospheric motion.

system (12)-(13). So, this might result in solutions that look “non-physical”. Our goal for the forthcoming studies is to investigate the more complex forms of the similarity solution like in (65) and we also aim to study numerically the complete nonlinear model (12)-(13) by using the perturbation analysis, in which the invariant solutions obtained in this paper will serve at the unperturbed solutions. We hope to find the perturbed solutions that could be associated with more physically relevant solutions like in schematic shown in Figure 9, which demonstrate the perturbation of atmospheric waves with various wave numbers.

Acknowledgements

We would like to thank Professor Derin Wysham (Department of Mathematics, Wenatchee Valley College) for his great help in attracting undergraduate students to this research project and useful discussions. Our grateful thanks are also extended to the Dean of STEM & Social Sciences at Wenatchee Valley College, Holly Bringman for providing us with all the necessary facilities for the research, great support, and encouragement.

NOTES

^{1}The symmetries of the system (14)-(15) can be reproduced from the results of the paper [38]. See also [42] and [43].

References

[1] 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

[2] 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

[3] 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

[4] 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

[5] 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

[6] Newman, P. (2010) Chemistry and Dynamics of the Antarctic Ozone Hole, in the Stratosphere: Dynamics, Transport, and Chemistry, Geophys. Monogr. Ser.

https://doi.org/10.1029/2009GM000873

[7] 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

[8] 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

[9] 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

[10] Baldwin, M. and Dunkerton, T. (2001) Stratospheric Harbingers of Anomalous Weather Regimes. Science, 294, 581-584.

https://doi.org/10.1126/science.1063315

[11] El Niño and La Niña (2016) New Zealand’s National Institute of Water and Atmospheric Research. Archived from the Original on 11 April 2016.

[12] 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.

[13] 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

[14] Sun, L., Wu, H. and Li, X. (2006) Our Understanding of Arctic Vortex. Chinese Journal of Polar Research, 18, 52-62.

[15] McCreary, J. (1976) Eastern Tropical Ocean Response to Changing Wind Systems: With Applications to El Niño. Journal of Physical Oceanography, 6, 632-645.

https://doi.org/10.1175/1520-0485(1976)006<0632:ETORTC>2.0.CO;2

[16] Murray, F. (1960) Dynamic Stability in the Stratosphere. Journal of Geophysical Research, 65, 3273-3305.

https://doi.org/10.1029/JZ065i010p03273

[17] 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

[18] Biebricher, A. (2015) Planetary Waves: The Biggest Coherent Pressure Structures in the Earth’s Atmosphere, Space for Science.

[19] 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

[20] 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

[21] Ibragimov, R. and Gunag, L. (2015) Splitting Phenomenon of a Higher-Order Shallow Water Theory Associated with a Longitudinal Planetary Wave. Dynamics of Atmospheres and Oceans, 69, 1-11.

https://doi.org/10.1016/j.dynatmoce.2014.10.003

[22] 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

[23] Schoeberl, M. and Hartmann, D. (1991) The Dynamics of the Stratospheric Polar Vortex and Its Relation to Springtime Ozone Depletions. Science, 251, 46-52.

https://doi.org/10.1126/science.251.4989.46

[24] Vassada, A., Horst, S., Kennedy, M. and Ingersoll, A. (2006) Cassini Imaging of Saturn: Southern Hemisphere Winds and Vorticities. Journal of Geophysical Research, 111, 5004-5017.

https://doi.org/10.1029/2005JE002563

[25] Stanford University (2018) Pulse of the Polar Vortex Revealed: A Key to Mapping Future Storms.

[26] 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

[27] Yang, X. (2016) Variability of the Northern Circumpolar Vortex and Its Association with Climate Anomaly in China. Acta Meteorologica Sinica, 26, 135-142.

[28] 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

[29] 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

[30] 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

[31] 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

[32] Lie, S. (1891) Vorlesungen uber Differentialgleichungen Mit Bekannten Infinitesimalen Transformationen, Bearbeited und herausgegeben von Dr. G. Scheffers, B.G. Teubner, Leizig.

[33] 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.

[34] Dickson, L.E. (1924) Differential Equations from the Group Standpoint. Annals of Mathematics, 25, 287-378.

https://doi.org/10.2307/1967773

[35] Ovsyannikov, L.V. (1962) Group Properties of Differential Equations. Siberian Branch, USSR Academy of Sciences, Novosibirsk.

[36] Ibragimov, N.H. (1983) Transformation Groups Applied to Mathematical Physics. Nauka, Moscow.

[37] 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

[38] 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

[39] Stephani, H. (1989) Differential Equations: Their Solution Using Symmetries. Cambridge Univ. Press, Cambridge.

https://doi.org/10.1017/CBO9780511599941

[40] Cantwell, B.J. (2002) Introduction to Symmetry Analysis. Cambridge University Press, Cambridge.

[41] Bluman, G.W. and Anco, S.C. (2002) Symmetry and Integration Methods for Differential Equations. Springer-Verlag, New York.

[42] Ibragimov, N. (2009) A Practical Course in Differential Equations and Mathematical Modelling. Higher Education Press, Beijing.

https://doi.org/10.1142/7573

[43] Kovalev, V., Pustovalov, V. and Senashov, S. (1993) Lie-Bäcklund Group for the Nonlinear Geometrical Optics Equations. Differential Equations, 29, 1751-1764.