The study of solitary and shallow water waves has gained a great significance to describe characteristics of nonlinear wave phenomena. This nonlinear feature arises due to interaction of propagating high and low amplitude waves. Due to highly nonlinear behavior, such wave phenomena are expressed by nonlinear partial differential equations (NPDEs). These NPDEs most often illustrate the dynamics of shallow water waves, interaction of propagating solitary waves, flow of heat and mass, ion acoustic waves in plasma physics and electromagnetic waves in optical fibers etc. [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41]. Now a days, numerous great efforts have been spent to investigate exact solutions of NPDEs and to understand the dynamics of associated nonlinear physical models. Thus, several adequate tools notably variable separation approach [4], extended mapping approach [6], rational expansion method [7], extended Fan sub-equation method [8], [9], Exp-function method [10], [18], $\left(\frac{G^{\prime }}{G}\right)$ -expansion method [11], Hirota bilinear method [12], [16], Bäcklund transformation method [14], [26], modified Kudryashov method [21], [22], [23], modified Sardar sub-equation [24], generalized auxiliary equation method [25] and Lie symmetry approach [5], [15], [27], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41] have been evolved for deep observations of these NPDEs.

To reveal the significant features, we are interested to enlighten some new invariant solutions of Broer-Kaup-Kupershmidt (BKK) equation

where $u$ and $v$ represent the height and horizontal velocity of water waves respectively. It represents the propagation of nonlinear and dispersive long gravity waves in horizontal directions $x$ and $y$ in shallow water of uniform depth. The BKK equation can be deduced from Kadomtsev-Petviashvili (KP) equation using inner parameter symmetry constraints [1]. While $x=y$ in Eq. (1.1) results in one dimensional Broer-Kaup equations, which describes dispersive long waves traveling in shallow water. The shallow water waves are the progressive waves, which classify marine environment, examine the ocean dynamics and model the equatorial tsunami waves. These waves arise due to action of wind on water and describe the long wavelength phenomena where horizontal length is greater than fluid depth. During propagation, these waves are affected by ocean bottom and thus orbital motion of water is disturbed. Therefore, it may cause underwater earthquakes, landslides, unimaginable destruction to the coastal ecology etc. Moreover, the internal waves arising from temperature variation can destruct marine vessels and can affect ecological as well as socio-economic environment.

A rich literature of BKK equation is reviewed here as Yu and Xin [2] analyzed PainlevȨ property and derived infinitely many symmetries of Eq. (1.1). Zhang et al. [3] used truncated PainlevȨ expansion method and obtained some soliton solutions of generalized BK equations. With the aid of variable separation approach, H.M. LI [4] found localized coherent structures of Eq. (1.1) and introduced its soliton solutions. H.C. Ma [5] derived Lie symmetries and then found exact dromion solutions. Thereafter, Fang and Zheng [6] studied BKK equation and derived nonelastic soliton solutions using extended mapping approach. Some soliton and periodic solutions were constructed by wang and Chen [7]. Zhang and Xia [8] proposed Fan sub-equation method and dummyTXdummy- used it to generated soliton and periodic solutions of Eq. (1.1). Continuing with same method, E. Yomba [9] investigated some general non-traveling wave solutions. Using Exp-function method via Riccati equation, some exact solutions of BKK equation were derived by S. Zhang [10]. Later on, Zhang et al. [11] employed $\left(\frac{G^{\prime }}{G}\right)$ -expansion method to Eq. (1.1) and listed some exact solutions. By means of symbolic computation and Hirota bilinear method, some soliton solutions were constructed by X.Y. Wen [12]. Moreover, Hu and Chen [13] derived nonlocal symmetries for BKK equation and studied its soliton and periodic solutions. By virtue of symbolic computation and Bell polynomial approach, Lan et al. [14] investigated some soliton solutions of Eq. (1.1). Furthermore, Kassem and Rashed [15] generated optimal system of Lie symmetries and derived some particular exact solutions. Some Wronskian solutions were introduced by Tang et al. [16] using Hirota method and bell polynomial approach. Recently, Kumar et al. used Lie symmetry method and found some exact solutions [17].

Being motivated with rich literature, we derive Lie symmetries and soliton solutions of BKK Eq. (1.1). Lie symmetry analysis under one parameter group of transformations ensures the invariance of system, which is used make the PDEs integrable and to reduce its order. It is also used to derive point symmetries of PDEs. Therefore, the given system of PDEs produces an equivalent system of ODEs and results in exact solutions. Based on characteristic of system and existing parameters, these solutions exhibit parabolic nature, soliton nature, intensive and non-intensive behavior of solitons. The solitons have wide applications in diverse physical contexts like plasma physics, quantum fluid dynamics and stratified fluid flows. Thus, Lie symmetry method serves as most applicable mathematical tool to find exact solutions having applications in various physical contexts.

The remaining texts are structured as: Section 2 presents the basic terms of symmetry method and invariance of BKK equation. The symmetry reductions and exact solutions are listed in Section 3. The derived results are analyzed in Section 4. The Section 5 is furnished with concluding remarks of findings in the end.

The section briefly defines the termilogy to exatract infinitesimal generators. The detailed methodology can be reviewed from literature [27], [28] and existing references. Let us consider the one-parameter ($ϵ$ ) transformations

where $x_1=t$ , $x_2=x$ , $x_3=y$ , $u_1=u$ , $u_2=v$ are variables, and ${\eta }^1=\tau$ , ${\eta }^2=\xi$ , ${\eta }^3=\eta$ , ${\psi }^1=\varphi$ , ${\psi }^2=\psi$ are infinitesimals designated to $t$ , $x$ , $y$ , $u$ , $v$ for maintaining the invariance of Eqs. ${\Delta }_1=0$ and ${\Delta }_2=0$ . The infinitesimal vector is $\begin{array}{cc}\hfill {\displaystyle }& \mathbf{v}=\tau {\partial }_t+\xi {\partial }_x+\eta {\partial }_y+\varphi {\partial }_u+\psi {\partial }_v.\hfill \end{array}$ Employing prolongation formulas $P{r}^{\left(3\right)}\mathbf{v}\left({\Delta }_1\right)=0$ and $P{r}^{\left(1\right)}\mathbf{v}\left({\Delta }_2\right)=0$ to Eq. (1.1), we have invariance condition as

where the extension coefficients ${\varphi }^x$ , ${\varphi }^y$ , ${\varphi }^{yt}$ , ${\varphi }^{xxy}$ , ${\psi }^t$ , ${\psi }^x$ , ${\psi }^{xx}$ , ${\psi }^{xy}$ are defined as $\begin{array}{ccc}\hfill {\varphi }^x& =& D_x\varphi -u_x{D}_x\xi -u_y{D}_x\eta -u_t{D}_x\tau ,\hfill \\ \hfill {\varphi }^y& =& D_y\varphi -u_x{D}_y\xi -u_y{D}_y\eta -u_t{D}_y\tau ,\hfill \\ \hfill {\psi }^t& =& D_t\psi -u_x{D}_t\xi -u_y{D}_t\eta -u_t{D}_t\tau ,\hfill \\ \hfill {\psi }^x& =& D_x\psi -u_x{D}_x\xi -u_y{D}_x\eta -u_t{D}_x\tau ,\hfill \\ \hfill {\psi }^{xx}& =& D_x^2\psi -u_x{D}_x^2\xi -u_y{D}_x^2\eta -u_t{D}_x^2\tau \hfill \\ & & -2u_{xx}{D}_x\xi -2u_{xy}{D}_x\eta -2u_{xt}{D}_x\tau ,\hfill \\ \hfill {\psi }^{xy}& =& D_x{D}_y\psi -u_x{D}_x{D}_y\xi -u_y{D}_x{D}_y\eta -u_t{D}_x{D}_y\tau \hfill \\ & & -u_{xy}{D}_x\xi -u_{xx}{D}_y\xi -u_{yy}{D}_x\eta -u_{xy}{D}_y\eta -u_{yt}{D}_x\tau -u_{xt}{D}_y\tau ,\hfill \\ \hfill {\varphi }^{yt}& =& D_y{D}_t\varphi -u_x{D}_y{D}_t\xi -u_y{D}_y{D}_t\eta -u_t{D}_y{D}_t\tau \hfill \\ & & -u_{xt}{D}_y\xi -u_{yt}{D}_y\eta -u_{tt}{D}_y\tau -u_{xy}{D}_t\xi -u_{yy}{D}_t\eta -u_{yt}{D}_t\tau ,\hfill \\ \hfill {\varphi }^{xxy}& =& D_x^2{D}_y\varphi -u_{xxx}{D}_y\xi -2u_{xxy}{D}_x\xi -2u_{xyy}{D}_x\eta -2u_{xyt}{D}_x\tau \hfill \\ & & -u_{xxy}{D}_y\eta -u_{xxt}{D}_y\tau -2u_{xx}{D}_x{D}_y\xi -2u_{xy}{D}_x{D}_y\eta -2u_{xt}{D}_x{D}_y\tau \hfill \\ & & -u_x{D}_x^2{D}_y\xi -u_{xy}{D}_x^2\xi -u_{yy}{D}_x^2\eta -u_{yt}{D}_x^2\tau -u_y{D}_x^2{D}_y\eta -u_t{D}_x^2{D}_y\tau ,\hfill \end{array}$ with total derivative operators $D_x$ , $D_y$ , $D_t$ .

The Eqs. (2.1) and (2.2) along with given expressions yield determining equations as $\begin{array}{ccc}& & {\xi }_u=0,{\xi }_v=0,{\xi }_y=0,{\eta }_x=0,{\eta }_t=0,{\eta }_u=0,{\eta }_v=0,\hfill \\ & & {\tau }_x=0,{\tau }_y=0,{\tau }_u=0,{\tau }_v=0,\hfill \\ & & 2{\xi }_x-{\tau }_t=0,2\varphi +u{\tau }_t-{\xi }_t=0,2\psi +v{\tau }_t+2v{\eta }_y=0.\hfill \end{array}$ This system generates the desired infinitesimals $\begin{array}{ccc}\hfill \xi & =& \frac{1}{2}x\overline{f_2}\left(t\right)+f_3\left(t\right),\eta =f_1\left(y\right),\tau =f_2\left(t\right),\hfill \\ \hfill \varphi & =& \frac{1}{4}(-2u\overline{f_2}\left(t\right)+x\overline{\overline{f_2}}\left(t\right)+2\overline{f_3}\left(t\right)),\hfill \\ \hfill \psi & =& -\frac{1}{2}v(2\overline{f_1}\left(y\right)+\overline{f_2}\left(t\right)),\hfill \end{array}$ with arbitrary functions $f_1\left(y\right)$ , $f_2\left(t\right)$ , $f_3\left(t\right)$ while the bar denotes derivative.

Hence, algebra of infinitesimal symmetries are generated by following vectors

$\begin{array}{cc}& {\mathbf{v}}_{\mathbf{1}}\left(f_1\right)=f_1{\partial }_y-v\overline{f_1}{\partial }_v,\hfill \\ & {\mathbf{v}}_{\mathbf{2}}\left(f_2\right)=\frac{x}{2}\overline{f_2}{\partial }_x+f_2{\partial }_t-\frac{u}{2}\overline{f_2}{\partial }_u+\frac{x}{4}\overline{\overline{f_2}}{\partial }_u-\frac{1}{2}v\overline{f_2}{\partial }_v,\hfill \\ & {\mathbf{v}}_{\mathbf{3}}\left(f_3\right)=f_3{\partial }_x+\frac{1}{2}\overline{f_3}{\partial }_u.\hfill \end{array}$ The commutative relations of above vectors are determined as

$\begin{array}{ccc}& & [{\mathbf{v}}_{\mathbf{1}},{\mathbf{v}}_{\mathbf{2}}]=0,[{\mathbf{v}}_{\mathbf{1}},{\mathbf{v}}_{\mathbf{3}}]=0,[{\mathbf{v}}_{\mathbf{2}},{\mathbf{v}}_{\mathbf{3}}]={\mathbf{v}}_{\mathbf{3}}(f_2\overline{f_3}-\frac{1}{2}\overline{f_2}{f}_3),\\ & & [{\mathbf{v}}_{\mathbf{2}},{\mathbf{v}}_{\mathbf{1}}]=0,[{\mathbf{v}}_{\mathbf{3}},{\mathbf{v}}_{\mathbf{1}}]=0,[{\mathbf{v}}_{\mathbf{3}},{\mathbf{v}}_{\mathbf{2}}]=-{\mathbf{v}}_{\mathbf{3}}(f_2\overline{f_3}-\frac{1}{2}\overline{f_2}{f}_3),\\ & & [{\mathbf{v}}_{\mathbf{1}}\left(g_1\right),{\mathbf{v}}_{\mathbf{1}}\left(h_1\right)]={\mathbf{v}}_{\mathbf{1}}(g_1{\bar{h}}_1-h_1\overline{g_1}),[{\mathbf{v}}_{\mathbf{2}}\left(g_2\right),{\mathbf{v}}_{\mathbf{2}}\left(h_2\right)]\\ & & ={\mathbf{v}}_{\mathbf{2}}(g_2{\bar{h}}_2-h_2{\bar{g}}_2),[{\mathbf{v}}_{\mathbf{3}}\left(g_3\right),{\mathbf{v}}_{\mathbf{3}}\left(h_3\right)]=0,\end{array}$ which constitutes infinite dimensional Lie algebra of symmetries.

To derive invariant solutions of BKK equation, the characteristic equation is furnished by

To continue integration further, the following two cases are taken under consideration as:

Case (I): If ${\displaystyle f_2\left(t\right)=\frac{1}{a}}$ , where $a\ne 0$ , then Eq. (3.1) reforms to
$\begin{array}{cc}& {\displaystyle \frac{dx}{f_3\left(t\right)}=\frac{dy}{f_1\left(y\right)}=\frac{adt}{1}=\frac{2du}{\overline{f_3}\left(t\right)}=\frac{-dv}{v\overline{f_1}\left(y\right)}.}\hfill \end{array}$ Then, group invariants of BKK equation yield

with symmetry variables $\begin{array}{cc}\hfill {\displaystyle r=}& {\displaystyle x-a\int f_3\left(t\right)dt\text{and}s=\int \frac{dy}{f_1\left(y\right)}-at.}\hfill \end{array}$ Thereupon, first symmetry reduction of BKK equation results

For further reduction, infinitesimal generators $\tilde{\xi }$ , $\tilde{\eta }$ , $\tilde{\varphi }$ and $\tilde{\psi }$ are

$\begin{array}{cc}\hfill {\displaystyle }& \tilde{\xi }=a_{1}r+a_2,\tilde{\eta }=2a_{1}s+a_3,\tilde{\varphi }=-a_{1}G,\tilde{\psi }=-3a_{1}H\hfill \end{array}$ with constants $a_1$ , $a_2$ and $a_3$ . So, one gets

Here, two intermediary cases are illustrated as:

Case (I$\mathbf{a}$ ): For $a_1\ne 0$ , Eq. (3.5) turns to $\begin{array}{ccc}\hfill {\displaystyle \frac{dr}{r+A_1}}& =& \frac{ds}{2(s+A_2)}=-\frac{dG}{G}=-\frac{dH}{3H}.\hfill \\ & & (A_1=\frac{a_2}{a_1}\text{and}{A}_2=\frac{a_3}{2a_1})\hfill \end{array}$ Eventually, the group invariants are $G=\frac{G_1\left(X\right)}{{(s+A_2)}^{\frac{1}{2}}}\text{and}H=\frac{H_1\left(X\right)}{{(r+A_1)}^3},$ with symmetry variable ${\displaystyle X=\frac{{(r+A_1)}^2}{s+A_2}}$ .

Substitution of above invariants into Eqs. (3.3) and (3.4) yields

The desired primitives are $G_1=-\frac{a}{2}{X}^{\frac{1}{2}}\text{and}{H}_1=c_1{X}^2$
where $c_1$ is constant.

Consequently, one derives the solutions

Case (Ib): For $a_1=0$ and $a_3\ne 0$ , Eq. (3.5) reforms to

$\begin{array}{cc}\hfill {\displaystyle }& \frac{dr}{A_3}=\frac{ds}{1}=\frac{dG}{0}=\frac{dH}{0}.(A_3=\frac{a_2}{a_3})\hfill \end{array}$ It provides

$G=G_2\left(Y\right)\text{and}H=H_2\left(Y\right)$ with variable ${\displaystyle Y=r-A_{3}s}$ .

Availing above invariants into Eq. (3.6) and (3.7), we have

The primitives are read as $\begin{array}{cc}\hfill {\displaystyle }& {\displaystyle G_2=-\frac{a{A}_3}{2}-\frac{1}{Y+c_3}\text{and}{H}_2=c_2(Y+c_3)}\hfill \end{array}$ which rises to

Case (II): For ${\displaystyle f_1\left(y\right)=b}$ and ${\displaystyle f_3\left(t\right)=c\overline{f_2}\left(t\right)}$ , Eq. (3.1) reforms to

$\begin{array}{ccc}\hfill {\displaystyle \frac{2dx}{(x+2c)\overline{f_2}\left(t\right)}}& =& \frac{dy}{b}=\frac{dt}{f_2\left(t\right)}=\frac{4du}{(x+2c)\overline{\overline{f_2}}\left(t\right)-2u\overline{f_2}\left(t\right)}\hfill \\ & =& \frac{-2dv}{v\overline{f_2}\left(t\right)}.\hfill \end{array}$ Then, group invariants yields

with symmetry variables $\begin{array}{c}\hfill {\displaystyle r=\frac{x+2c}{f_2^{\frac{1}{2}}\left(t\right)}\text{and}s=y-b\int \frac{1}{f_2\left(t\right)}dt.}\end{array}$ Thus, first symmetry reduction produces

To proceed further, infinitesimal generators $\tilde{\xi }$ , $\tilde{\eta }$ , $\tilde{\varphi }$ and $\tilde{\psi }$ are $\begin{array}{cc}\hfill {\displaystyle }& \tilde{\xi }=a_{4}r+a_5,\tilde{\eta }=2a_{4}s+a_6,\tilde{\varphi }=-a_{4}G,\tilde{\psi }=-3a_{4}H,\hfill \end{array}$
where $a_4$ , $a_5$ and $a_6$ are constants.

Therefore, characteristic equation is

Again, two intermediary cases are illustrated as:

Case (IIa): For $a_4\ne 0$ , Eq. (3.17) generates $\begin{array}{cc}\hfill {\displaystyle }& \frac{dr}{r+A_4}=\frac{ds}{2(s+A_5)}=-\frac{dG}{G}=-\frac{dH}{3H}\hfill \end{array}$
where ${\displaystyle A_4=\frac{a_5}{a_4}}$ and ${\displaystyle A_5=\frac{a_6}{2a_4}}$ . Then, group invariants are obtained as

$G=\frac{G_1\left(X\right)}{{(s+A_5)}^{\frac{1}{2}}}\text{and}H=\frac{H_1\left(X\right)}{{(r+A_4)}^3}$ with symmetry variable ${\displaystyle X=\frac{{(r+A_4)}^2}{s+A_5}}$ .

Using above values into system (3.15) and (3.16), one gets

The primitives of these are $G_1=-\frac{b}{2}{X}^{\frac{1}{2}}\text{and}{H}_1=c_4{X}^2.\left(c_4\text{being}\text{arbitrary}\text{constant.}\right)$

Consequently, one derives the solution

Case (IIb): For $a_4=0$ and $a_6\ne 0$ , Eq. (3.17) transforms to $\begin{array}{c}\hfill {\displaystyle \frac{dr}{A_6}=\frac{ds}{1}=\frac{dG}{0}=\frac{dH}{0}.(A_6=\frac{a_5}{a_6})}\end{array}$ It leads to

$G=G_2\left(Y\right)\text{and}H=H_2\left(Y\right)$ with symmetry variable ${\displaystyle Y=r-A_{6}s}$ .

Availing these invariants into Eqs. (3.18) and (3.19), we have

It has primitives

$G_2=-\frac{b{A}_6}{2}-\frac{1}{Y+c_6}\text{and}{H}_2=c_5(Y+c_6)$ and lead to

where ${\displaystyle c_5}$ and ${\displaystyle c_6}$ are constants.

A comparison of derived solutions (3.8), (3.9), (3.12), (3.13), (3.20), (3.21), (3.24) and (3.25) with previous results is presented here. The derived solutions involve the arbitrary functions $f_1\left(y\right)$ , $f_2\left(t\right)$ , $f_3\left(t\right)$ and arbitrary constants as well. The derived solutions are entirely new and general than previous finding [15], [17], [19] as

● Kassem and Rashed [15] derived their solutions by restricting the arbitrary functions and taking their values as constants while the solutions listed in this work contain several arbitrary functions.

● Kumar et al. [17] derived symmetry reductions and exact solutions for the single vector or combination of two vector fields only, while the present work deals with the general vector field $\mathbf{v}={\mathbf{v}}_{\mathbf{1}}\left(f_1\right)+{\mathbf{v}}_{\mathbf{2}}\left(f_2\right)+{\mathbf{v}}_{\mathbf{3}}\left(f_3\right)$ . Therefore, the derived solutions involves more arbitrary functions existed in infinitesimals, which may be helpful to derives several physical phenomena and enhance the significance of results.

● Li et al. [19] derived the results involving arbitrary function of $y$ only, while the solutions listed in the present work contain arbitrary function of $y$ and $t$ . Moreover, the present results are entirely different from Li et al. [19].

Thus, the comparison with previous results ensures the novelty of results and existence of arbitrary functions shows enrichment of present work.

The graphical interpretation of invariant solutions encourages to reveal significant features physical phenomena. It demonstrates the precise informations to predict phenomena dynamics. This segment reveals the physical behavior of the phenomena associated with invariant solutions (3.8), (3.9), (3.12), (3.13), (3.20), (3.21), (3.24) and (3.25) using numerical simulation. The solutions involve arbitrary functions $f_1\left(y\right)$ , $f_2\left(t\right)$ , $f_3\left(t\right)$ and several arbitrary constants so simulation is performed for their appropriate values $f_1\left(y\right)=\frac{1}{2(\alpha +y)}$ , $\alpha =0.8351$ , $a=0.5825$ , $b=0.8929$ , $c=0.7032$ , $A_1=0.5827$ , $A_2=0.8549$ , $A_3=0.587$ , $A_4=0.2176$ , $A_5=0.251$ and rest values are listed in adjacent figures. Due to existence of these functions and constants, the dynamical system has rich physical structures. Consequently, doubly soliton, multi-soliton, parabolic nature and soliton fission behavior have been analyzed. Solitons are nonlinear solitary waves, which generally preserve their identities after mutual collision and show elastic behavior. But for instance, interactions of solitons may not be elastic and result into soliton fusion or soliton fission. The detailed analysis of solution profiles is performed as follows:

Fig. 1: The surface plots of solution profiles $u$ and $v$ for Eqs. (3.8) and (3.9) are shown in this figure. The components $u$ and $v$ show doubly and multisolition profiles for $f_3\left(t\right)=2(\beta +t)$ , $\beta =0.8954$ , $c_1=0.0348$ at $t=0.7551$ . The interaction of solitons is observed completely elastic with high peaks and deep tails under the plane wave.

Fig. 2: The spatio-temporal profiles for wave components $u$ and $v$ listed in Eqs. (3.12) and (3.13) are traced with adjacent figure. The values of functions and parameters are considered as $f_3\left(t\right)=2(\beta +t)$ , $\beta =0.8954$ , $c_2=0.9421$ , $c_3=0.9564$ at $y=0.1319$ . The component $u$ shows soliton profile while $v$ shows parabolic nature. Thus during propagation, the waves interact and result into nonlinear parabolic nature.

Figs. 3 and 4: The Fig. 3 reveals the propagation of multisoliton and doubly soliton for Eqs. (3.20) and (3.21). The appropriate values of functions and constants are taken $f_2\left(t\right)=(\gamma +t)$ , $\gamma =0.5557$ , $c_4=0.1844$ . The interaction of multisoliton and doubly soliton is observed at different values of $t$ and shown at $t=0.3174$ in this figure, which results into multisoliton after mutual collision. The $y-u$ and $y-v$ views show that waves attain highest amplitude around origin. Moreover, treating a periodic function $f_2\left(t\right)=\sec (\delta +t)$ with $\delta =0.0319$ , solution profile preserves the similar soliton nature shown in the Fig. 4. However, the amplitude of waves vary with respect to time, different values of functions and parameters as well.

Fig. 5: The spatio-temporal profile for Eqs. (3.20) and (3.21) is traced via this figure. It reveals soliton fission phenomena for appropriate values $f_2\left(t\right)=(\gamma +t)$ , $\gamma =0.5557$ , $c_4=0.1844$ at $y=0.9597$ . The fission of one soliton into two is exhibited in $t-u$ view and shows new phenomena. Thus, it is pointed out that physical phenomena preserve its behavior with respect to spatial variables while the behavior is changed for space-time variables.

Fig. 6: The wave components expressed in Eqs. (3.24) and (3.25) are analyzed through this spatio-temporal surface plot, which reflects the elastic behavior multisoliton at $y=0.7922$ . The values to the constants are provided as $f_2\left(t\right)=\sec (\delta +t)$ , $\delta =0.0319$ , $c_5=0.9594$ , $c_6=0.6557$ . The waves with different amplitudes interact during propagation and ultimately re-merge to its own shape, which is represented in $t-u$ views.

In this work, we study the BKK equation and analyze its physical significance in the shallow water of uniform depth. The BKK equation illustrates the nonlinear and dispersive long gravity waves propagating in two horizontal directions. The invariance of system under one parameter transformations leads to infinitesimal generators. The commutative relations of infinitesimal vectors are derived, which constitute infinite dimensional Lie algebra of symmetries. Thereafter, Lie symmetry method is employed to derive symmetry reductions and generate exact solutions listed in Eqs. (3.8), (3.9), (3.12), (3.13), (3.20), (3.21), (3.24) and (3.25). These solutions retain arbitrary functions involved in infinitesimals as well as some free parameters. Due the existence of arbitrary functions and parameters, these finding have rich physical structures and can be used to analyze the dynamics of real phenomena. The physical structures of solutions are analyzed using numerical simulation, which show the nonlinear parabolic nature and characteristics of soliton interactions. Thus, Lie symmetry method can be treated as an effective tool to handle system of nonlinear PDEs.

Authors declare that they have no conflict of interest.