3.1.1. Decomposition of collapse sub-processes

An initial profile of bulge is depicted by an orange line in figure 5(a). As time proceeds, the whole of the initial bulge is collapsing and spreading outwards due to gravity, with the collapse occurring from its axisymmetric centre towards the surrounding stationary film. Hence, this process is named the ‘collapse process’. During the collapse process, it is generally expected that the hydrostatic pressure at the centre of the bulge is greater than that of the surrounding zone. Under the influence of gravity, the bulge collapses and ‘squeezes’ outward. Obviously, after a certain period, as depicted in figure 5(a), the centre bulge height, which corresponds to the peak of the initial bulge (orange line), descends below the surrounding zone, forming what can be referred to as a ‘basin’ (green line). This phenomenon indicates that fluids move from the centre to the surrounding zone, resulting in a new wavefront around it. To facilitate analyses and discussions, we designate the shape where the bulge’s centre height is slightly below the height of the surrounding zone as the ‘transition shape’. It is noticeable that the process from the start of collapse to the formation of the ‘transition shape’ is influenced by the surrounding flow fluid and gradually squeezes a new wavefront.

Figure 5.

Collapse process. (a) Sub-process I: squeezing process. (b) Sub-process II: diffusion process.

Once the centre height of the bulge falls below the surrounding zone, from the perspective of energy, the gravitational potential energy at the centre of the bulge is converted into the kinetic energy of the newly extruded wavefront. This wavefront, acting as a leading wave, propagates radially outwards. It is generally considered that the new wavefront corresponds to propagation with a high Re, while the diffusion behaviour of the centre has a low Re of the order of $O(1)$ (Kakutani & Matsuuchi Reference Kakutani and Matsuuchi1975). Therefore, the newly formed wavefront will continue to move with a significantly higher overall velocity, as shown in figure 5(b). At the stage from the formation of ‘transition shape’ to the point where the centre height of the bulge reaches the stationary water level, the centre of the bulge diffuses out without being influenced by the surrounding propagating leading wave. For further analysis, we will divide the collapse process into two distinct sub-processes. Sub-process I begins with the initial shape of the bulge and ends with the ‘transition shape’, during which the decrease in central wave height is significantly influenced by the surrounding fluid and squeezes out new wavefronts. Therefore, this is referred to as the ‘squeezing process’. While sub-process II commences with the ‘transition shape’ and extends until the end of the collapse process, where the centre of the bulge diffuses without influence from the surrounding new wavefronts, thus constituting a ‘diffusion process’.

To rigorously justify that the propagating leading wave will move with a significantly higher overall velocity and will not affect the diffusion process, we monitor the overall velocities of these two processes as shown in figure 6. In the diffusion process at time instance $t_1$ in figure 6(a), we use the wave height $H'$ at the symmetrical axis to reconstruct a corresponding solitary wave with the identical height. Then, we find the measurement point on the solitary wave with a half-wave height. We use the radius velocity component of the measurement point to represent the overall velocity of the diffusion process $U_d$ . As for the propagation process in figure 6(b), its overall velocity can be easily evaluated by the phase velocity $U_p = ({r_3-r_2})/({t_3-t_2})$ . Our DNS results indicate that the overall velocities of these two processes have a relationship $U_d/U_p\sim O(10^{-2})$ , which prove that the diffusion process is an independent process.

Figure 6.

Overall velocity measurements. (a) Diffusion process. (b) Propagation process.

3.1.2. Scaling and perturbation analyses

Due to the interaction of the surrounding fluid, the squeezing process is a relatively complex phenomenon, making it difficult to characterise the wave height variation using a simple equation. In contrast, the diffusion at the centre of the bulge is not significantly influenced by the surrounding new wavefront during the diffusion process. Therefore, we will conduct a detailed modelling analysis specifically for the diffusion process. In the diffusion process, our main objective is to derive a simplified equation for the surface elevation under the assumption of weak nonlinearity and long wavelength. Similar to the approach used by Kakutani & Matsuuchi (Reference Kakutani and Matsuuchi1975) for studying a 2-D diffusion problem, the slow radial variable and slow time variable are adopted: $\tilde \xi = \varepsilon \tilde r;\tilde \tau = {\varepsilon ^2}\tilde t$ , which are substituted into (3.1) and (3.2). These transformations are required to ensure that $\alpha$ is of order $\varepsilon$ . Here, $\varepsilon = \eta /h$ represents the amplitude parameter, measuring nonlinearity, and $\alpha = h/\lambda$ represents the long-wave parameter, measuring dispersion. The undisturbed steady state without the bulge collapse is $\tilde u = \tilde w = 0$ and $\tilde p = {\tilde p_\infty } - \tilde z$ . By introducing (3.3), the form of asymptotic multiple-scale expansions, into the governing equations and combining with boundary conditions, a sequence of equations and boundary conditions for each order of $\varepsilon$ can be obtained. It is noted that any field quantity $\tilde \phi (\tilde \xi ,\tilde z = \tilde \eta ,\tilde \tau )$ at the free surface boundary condition can be expressed as a Taylor expansion equation (3.4) at the position $\tilde z_0$ to match the multiple-scale expansions in the governing equations.

(3.3) \begin{align}& \left \{ \begin{array}{l} \tilde u\big( {\tilde \xi ,\tilde z,\tilde \tau } \big) = \sum \limits _{n = 0}^\infty {{\varepsilon ^n}{{\tilde u}_n}\big( {\tilde \xi ,\tilde z,\tilde \tau } \big)} ;\quad \tilde w\big( {\tilde \xi ,\tilde z,\tilde \tau } \big) = \sum \limits _{n = 0}^\infty {{\varepsilon ^n}{{\tilde w}_n}\big( {\tilde \xi ,\tilde z,\tilde \tau } \big)}, \\ \tilde p\big( {\tilde \xi ,\tilde z,\tilde \tau } \big) = {{\tilde p}_\infty } - \tilde z + \sum \limits _{n = 0}^\infty {{\varepsilon ^n}{{\tilde p}_n}\big( {\tilde \xi ,\tilde z,\tilde \tau } \big)} ;\quad \tilde \eta \big( {\tilde \xi ,\tilde \tau } \big) = \sum \limits _{n = 0}^\infty {{\varepsilon ^n}{{\tilde \eta }_n}\big( {\tilde \xi ,\tilde \tau } \big)}, \end{array} \right . \end{align}

(3.4) \begin{align}&\qquad\qquad \tilde \phi \big( {\tilde \xi ,\tilde z = \tilde \eta ,\tilde \tau } \big) = { {\tilde \phi } \Bigg|_{\tilde z = {{\tilde z }_0}}} + {\left . {\frac {{\partial \tilde \phi }}{{\partial \tilde z}}} \right |_{\tilde z = {{\tilde z }_0}}}\left ( {\sum \limits _{n = 1}^\infty {{\varepsilon ^n}{{\tilde \eta }_n}\big( {\tilde \xi ,\tilde \tau } \big)} } \right ) + \ldots \\[12pt] \nonumber \end{align}

Due to the prevailing effect of viscous dissipation, the weaker nonlinearity cannot balance the dominant viscous dissipation, the position of Taylor expansion is selected as $\tilde z_0 = {\tilde \eta }_0$ for diffusion process in this section and $\tilde z_0 = 0$ for propagation process in next section (Kakutani & Matsuuchi Reference Kakutani and Matsuuchi1975). After substituting these expressions (3.3) and (3.4), the governing equation takes the form:

(3.5) \begin{equation} \left \{ \begin{array}{l} \displaystyle {\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + \varepsilon \left(\frac {{{{\tilde u}_0}}}{{\tilde \xi }} + \frac {{\partial {{\tilde u}_0}}}{{\partial \tilde \xi }} + \frac {{\partial {{\tilde w}_1}}}{{\partial \tilde z}}\right) + {\varepsilon ^2}\left(\frac {{{{\tilde u}_1}}}{{\tilde \xi }} + \frac {{\partial {{\tilde u}_1}}}{{\partial \tilde \xi }} + \frac {{\partial {{\tilde w}_2}}}{{\partial \tilde z}}\right) = 0} \\[13pt] \displaystyle {{{\tilde w}_0}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}} - \frac {1}{{{{Re}} }}\frac {{{\partial ^2}{{\tilde u}_0}}}{{\partial {{\tilde z}^2}}} + \varepsilon \left({{\tilde u}_0}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde \xi }} + {{\tilde w}_0}\frac {{\partial {{\tilde u}_1}}}{{\partial \tilde z}} + {{\tilde w}_1}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}} + \frac {{\partial {{\tilde p}_0}}}{{\partial \tilde \xi }} - \frac {1}{{{{Re}} }}\frac {{{\partial ^2}{{\tilde u}_1}}}{{\partial {{\tilde z}^2}}}\right)} \\[13pt] \displaystyle { + {\varepsilon ^2}\left(\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde \tau }} + {{\tilde u}_0}\frac {{\partial {{\tilde u}_1}}}{{\partial \tilde \xi }} + {{\tilde u}_1}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde \xi }} + {{\tilde w}_0}\frac {{\partial {{\tilde u}_2}}}{{\partial \tilde z}} + {{\tilde w}_1}\frac {{\partial {{\tilde u}_1}}}{{\partial \tilde z}} + {{\tilde w}_2}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}} + \frac {{\partial {{\tilde p}_1}}}{{\partial \tilde \xi }}\right .} \\[13pt] \displaystyle { \left . - \frac {1}{{{{Re}} }}\frac {{{\partial ^2}{{\tilde u}_2}}}{{\partial {{\tilde z}^2}}} + \frac {1}{{{{Re}} }}\frac {{{{\tilde u}_0}}}{{{{\tilde \xi }^2}}} - \frac {1}{{{{Re}} }}\frac {1}{{\tilde \xi }}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde \xi }} + \frac {1}{{{{Re}} }}\frac {{{\partial ^2}{{\tilde u}_0}}}{{\partial {{\tilde \xi }^2}}}\right) = 0} \\[13pt] \displaystyle {{{\tilde w}_0}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + \frac {{\partial {{\tilde p}_0}}}{{\partial \tilde z}} - \frac {1}{{{{Re}} }}\frac {{{\partial ^2}{{\tilde w}_0}}}{{\partial {{\tilde z}^2}}} + \varepsilon \left({{\tilde u}_0}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde \xi }} + {{\tilde w}_1}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + {{\tilde w}_0}\frac {{\partial {{\tilde w}_1}}}{{\partial \tilde z}} + \frac {{\partial {{\tilde p}_1}}}{{\partial \tilde z}} - \frac {1}{{{{Re}} }}\frac {{{\partial ^2}{{\tilde w}_1}}}{{\partial {{\tilde z}^2}}}\right)} \\[13pt] \displaystyle {+ {\varepsilon ^2}\left(\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde \tau }} + {{\tilde u}_1}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde \xi }} + {{\tilde u}_0}\frac {{\partial {{\tilde w}_1}}}{{\partial \tilde \xi }} + {{\tilde w}_1}\frac {{\partial {{\tilde w}_1}}}{{\partial \tilde z}} + {{\tilde w}_0}\frac {{\partial {{\tilde w}_2}}}{{\partial \tilde z}} + {{\tilde w}_2}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}}\right .}\\[13pt] \displaystyle {\left . + \frac {{\partial {{\tilde p}_2}}}{{\partial \tilde z}} - \frac {1}{{{{Re}} }}\frac {{{\partial ^2}{{\tilde w}_2}}}{{\partial {{\tilde z}^2}}} - \frac {1}{{{{Re}} }}\frac {1}{{\tilde \xi }}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde \xi }} - \frac {1}{{{{Re}} }}\frac {{{\partial ^2}{{\tilde w}_0}}}{{\partial {{\tilde \xi }^2}}}\right) = 0}. \end{array} \right . \end{equation}

Correspondingly, the boundary condition becomes

(3.6) \begin{equation} \left \{ \begin{array}{l} \displaystyle {\tilde u = \tilde w = 0\quad at\quad \tilde z = 0} \\[5pt] \displaystyle {{{\tilde w}_0} + {\varepsilon ^1}\left({{\tilde \eta }_1}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + {{\tilde w}_1} - {{\tilde u}_0}\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }}\right) + {\varepsilon ^2}\left[{{\tilde \eta }_2}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + {{\tilde \eta }_1}\frac {{\partial {{\tilde w}_1}}}{{\partial \tilde z}} + {{\tilde w}_2}\right .} \\[10pt] \displaystyle {\left .-\, \frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \tau }} - {{\tilde u}_0}\frac {{\partial {{\tilde \eta }_1}}}{{\partial \tilde \xi }} - \left({{\tilde \eta }_1}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}} + {{\tilde u}_1}\right)\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }}\right] = 0\quad at\quad \tilde z = 1 + {{\tilde \eta }_0}} \\[10pt] \displaystyle \frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}} + {{\varepsilon}^1} \left[\frac {{\partial \left({{\tilde \eta }_1}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}} + {{\tilde u}_1}\right)}}{{\partial \tilde z}} + \frac {{\partial {{\tilde w}_0}}}{{\partial \tilde \xi }}+ 2\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}}\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }}\right] \\[10pt] \displaystyle { +\, {\varepsilon ^2}\left\{ \frac {{\partial \left({{\tilde \eta }_2}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}} + {{\tilde \eta }_1}\frac {{\partial {{\tilde u}_1}}}{{\partial \tilde z}} + {{\tilde u}_2}\right)}}{{\partial \tilde z}} + \frac {{\partial \left({{\tilde \eta }_1}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + {{\tilde w}_1}\right)}}{{\partial \tilde \xi }} - \frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}}{\left(\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }}\right)^2}\right .} \\[15pt] \displaystyle {\left . +\, 2 \left[\frac {{\partial \left({{\tilde \eta }_1}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + {{\tilde w}_1}\right)}}{{\partial \tilde z}} - \frac {{\partial {{\tilde u}_0}}}{{\partial \tilde \xi }}\right]\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }} + 2\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}}\frac {{\partial {{\tilde \eta }_1}}}{{\partial \tilde \xi }}\right\} = 0\quad at\quad \tilde z = 1 + {{\tilde \eta }_0}} \\[15pt] \displaystyle { -\, {{\tilde p}_0} + {{\tilde \eta }_0} + \frac {2}{{{{Re}} }}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + {\varepsilon ^1}\left[\frac {2}{{{{Re}} }}\frac {{\partial \left({{\tilde \eta }_1}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + {{\tilde w}_1}\right)}}{{\partial \tilde z}} - \left({{\tilde \eta }_1}\frac {{\partial {{\tilde p}_0}}}{{\partial \tilde z}} + {{\tilde p}_1} - {{\tilde \eta }_1}\right)\right .} \\[20pt] \displaystyle {\left . -\, \frac {2}{{{{Re}} }}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}}\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }}\right] + {\varepsilon ^2}\left\{ {\left(\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }}\right)^2}\left ( { - {{\tilde p}_0} + {{\tilde \eta }_0}} \right ) - \left({{\tilde \eta }_2}\frac {{\partial {{\tilde p}_0}}}{{\partial \tilde z}} + {{\tilde \eta }_1}\frac {{\partial {{\tilde p}_1}}}{{\partial \tilde z}} + {{\tilde p}_2} - {{\tilde \eta }_2}\right) \right .} \\[13pt] \displaystyle {\left .+\, \frac {2}{{{{Re}} }}\frac {{\partial \left({{\tilde \eta }_2}\frac {{\partial {{\tilde w}_0}}}{{\partial \tilde z}} + {{\tilde \eta }_1}\frac {{\partial {{\tilde w}_1}}}{{\partial \tilde z}} + {{\tilde w}_2}\right)}}{{\partial \tilde z}} \right .} \\[13pt] \displaystyle {\left .-\, \frac {2}{{{{Re}} }}\left[\frac {{\partial \left({{\tilde \eta }_1}\frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}} + {{\tilde u}_1}\right)}}{{\partial \tilde z}}\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }} + \frac {{\partial {{\tilde u}_0}}}{{\partial \tilde z}}\frac {{\partial {{\tilde \eta }_1}}}{{\partial \tilde \xi }} + \frac {{\partial {{\tilde w}_0}}}{{\partial \tilde \xi }}\frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }}\right]\right\} = 0\quad at\quad \tilde z = 1 + {{\tilde \eta }_0}}. \end{array} \right . \end{equation}

With these multiple-scale expansion equations, the leading order solution of $\tilde \eta$ satisfying the following nonlinear diffusion equation is given as

(3.7) \begin{equation} \frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \tau }} = \frac {{{{Re}} }}{{3\tilde \xi }}\frac {\partial }{{\partial \tilde \xi }}\left [ {{{\left ( {1 + {{\tilde \eta }_0}} \right )}^3}\tilde \xi \frac {{\partial {{\tilde \eta }_0}}}{{\partial \tilde \xi }}} \right ] .\end{equation}

This equation is consistent with the diffusion equation derived by Nakaya (Reference Nakaya1974) in his study of the spread of a viscous droplet under gravity on a plane, which indicates that both belong to the same diffusion process. Furthermore, it also proves that during the diffusion process the surrounding leading wave exerts negligible influence on the centre of the bulge.

Figure 7.

Validations for diffusion process. (a) The DNS evolution of the wave height during the whole collapse process. All data points have been translated to start from an initial value of zero. The subscript $i$ denotes the selected initial value. Here, $h=0.3$ and $H=0.03$ are given with an initial bulge profile approximating a solitary wave. (b) Verification of scaling laws between wave height and time of diffusion process. All data points have been translated to start from an initial value of zero. The subscript $i$ denotes the selected initial value. The numbers following $h$ , $H$ , $W$ represent their values, expressed in SI units. The black solid lines, each with a slope of 45 degrees, indicate the trend given by the theoretically derived scaling law.

3.1.3. Normalisation of the scaling law

According to Landau & Lifshitz (Reference Landau and Lifshitz1987), (3.7) is solved by introducing the volume flow of bulge, $\tilde {q}$ , and the scaling law for the wave height at the central axis of the bulge over time is obtained:

(3.8) \begin{equation} {{\tilde \eta }_0} \sim {{\tilde q}^{1/4}}{{\tilde \tau }^{ - 1/4}}. \end{equation}

The evolution of wave height throughout the entire collapse process obtained through numerical simulation are depicted in figure 7(a). This figure distinctly illustrates that the height at the symmetrical axis deviates from the diffusion equation at the early stage. However, upon traversing the transition zone, the height attenuation satisfying the scaling law of (3.8), conforming its diffusion process. The marked transition zone in figure 7(a) corresponds to the transition shape depicted in figure 5(b). It proves that the entire collapse process comprises two distinct sub-processes: the squeezing process and the diffusion process. In the squeezing process, the bulge collapse is influenced by the ‘squeezed out’ new wavefront. Subsequently, in the diffusion process, the bulge diffuses freely.

Additionally, the diffusion process of bulge evolution with different $h$ and $H$ cases are conducted by DNS. These cases all fall within an initial wave profile that resembles a solitary wave. The results pertaining to the relationship between dimensionless wave height and time are depicted in figure 7(b). However, it is crucial to note that the initial volume flow during the diffusion process of the collapse process cannot be clearly defined, as it undergoes a squeezing stage. As a result, it cannot be normalised across different cases through initial volume flow. Since the dimensionless parameter $\varepsilon$ represents the relationship between initial wave height and water depth, which directly reflect the value of volume flow in the diffusion process, it is evident that the normalisation method for a bulge with an initial solitary wave profile is found to be associated with the parameter ${\varepsilon }^{- 5/4}$ . The obtained results exhibit a high degree of consistency with the solutions derived through theoretical analysis. It is worth noting that for the collapse of waves with different initial profiles, such as the collapse of an initial cosine wave, the scaling law in (3.8) is still satisfied, but the exponent of ‘ $\varepsilon$ ’ differs.

Figure 8.

Bulge-induced pressure model. (a) Comparison of wave shape between the solitary solution of the KdV equation and the solution of the propagation process of bulge evolution. (b) Comparison of wave decay between the solitary solution of the C-KdV equation and the solution of the propagation process of bulge evolution ( $h=0.1,H=0.02$ ). (c) Undular bore models in the 2-D axisymmetric condition.

3.2.1. Bulge-induced pressure model

In this section, we consider the propagation process. The new well-defined propagating wavefront continues to propagate along the radial direction. From the waveform obtained through numerical results, it can be observed that this new wavefront, also known as the leading wave, is followed by a series of nonlinear wave trains. As shown in figure 8(a), the solid black line represents the wave shape during the propagation process of bulge evolution. It is instructive to compare this wave shape with that of a solitary wave (dashed red line). An obvious nonlinear wave train is observed at the rear of the leading wave. Upon reviewing the waveform characteristics of the C-KdV equation (Miles Reference Miles1977), it is apparent that the propagation process in bulge evolution is very similar. However, the decay of the leading wave is not consistent. The wave height of the leading wave in the bulge evolution process decays faster compared with that in the C-KdV equation, as depicted in figure 8(b). One can infer from this difference that there is energy transfer from the leading wave to the wave train, then the leading wave gradually decays with time, whereas a solitary wave will propagate with a constant wave height and the C-KdV one will decay only because of a radial variation in the absence of the friction dissipation. Recent studies have introduced a model to visualise the effect of wind power on wave shape by incorporating a wind-induced surface pressure (Zdyrski & Feddersen Reference Zdyrski and Feddersen2021). Inspired by this work, a bulge-induced surface pressure $\displaystyle {{p_b} = {P_b}({1}/{r})({{(\partial r\eta)}}/{{(\partial r)}}})$ is explicitly introduced to represent the energy transfer from the leading wave to the wave train and added to the momentum equation along the radius direction as (3.9).

(3.9) \begin{equation} \frac {{\partial u}}{{\partial t}} + u\frac {{\partial u}}{{\partial r}} + w\frac {{\partial u}}{{\partial z}} = - \frac {1}{\rho }\frac {{\partial p}}{{\partial r}} + \frac {\mu }{\rho }\left ( { - \frac {u}{{{r^2}}} + \frac {1}{r}\frac {\partial }{{\partial r}}(r\frac {{\partial u}}{{\partial r}}) + \frac {{{\partial ^2}u}}{{\partial {z^2}}}} \right )\underbrace { - \frac {1}{\rho }\frac {{\partial {p_b}}}{{\partial r}}}_{{\mathrm{bulge\hbox{-}induced\thinspace pressure}}}. \end{equation}

The term ‘ $P_b$ ’ in bulge-induced surface pressure is assumed to be constructed using the well-known undular bore model, since there is a height difference between the leading wave and the wave train during the propagation process. It should be noted that there also exists a height difference between the upstream and downstream fluid in the undular bore model, as shown in figure 8(c). It is expected that based on the conservation laws of mass and momentum of the undular bore system, the rate of energy loss per unit volume per unit time through the bore can be evaluated. Lighthill (Reference Lighthill and Benjaminm1954) studied the energy dissipation rate in a 2-D condition and we extend this analysis to a 2-D axisymmetric condition. The conservation of mass and momentum are given as

(3.10) \begin{equation} \left \{ \begin{array}{l} {\mathrm{Mass\,conservation:}}\quad {u_1}{h_1}2\pi {r_1} = {u_2}{h_2}2\pi {r_2} = Q,\\ \displaystyle {{\mathrm{Momentum\, conservation:}}\quad \rho Q\left ( {{u_2} - {u_1}} \right ) = \frac {1}{2}\rho g{h_1} \cdot {h_1}2\pi {r_1} - \frac {1}{2}\rho g{h_2} \cdot {h_2}2\pi {r_2}}. \end{array} \right . \end{equation}

Hence, the rate of energy loss per unit volume per unit time can be expressed as

(3.11) \begin{equation} \displaystyle {\Delta E = \frac {1}{2}\rho \big( {u_1^2 - u_2^2} \big) + \rho g\left ( {{h_1} - {h_2}} \right )= \rho g\frac {{h_2^3\dfrac {{{r_2}}}{{{r_1}}} - 3h_2^2{h_1} + 3h_1^2{h_2} - h_1^3\dfrac {{{r_1}}}{{{r_2}}}}}{{4{h_1}{h_2}}}}. \end{equation}

Considering ${h_1} = h + {\eta _1}, {h_2} = h + {\eta _2}$ and a conserved quantity $\eta _1^{3/2}{r_1} = const = \eta _2^{3/2}{r_2}$ (Lamb Reference Lamb1932), then we can derive the transferring energy rate as

(3.12) \begin{equation} \displaystyle {\Delta E = \rho gh\frac {{{{\left ( {\dfrac {{{\eta _1}}}{{{\eta _2}}}} \right )}^{3/2}} - {{\left ( {\dfrac {{{\eta _1}}}{{{\eta _2}}}} \right )}^{ - 3/2}} + 3\left ( {\dfrac {{{\eta _1}}}{h} - \dfrac {{{\eta _2}}}{h}} \right )}}{{4\left ( {1 + \dfrac {{{\eta _1}}}{h}} \right )\left ( {1 + \dfrac {{{\eta _2}}}{h}} \right )}}}. \end{equation}

Due to $O({\eta _2}) \sim O({\eta _1}), ( {1 + {{{\eta _1}}}/{h}} ) \sim O(1)$ and $ ( {1 + {{{\eta _2}}}/{h}} ) \sim O(1)$ under the assumption of weak nonlinearity, $\Delta E$ is approximately ${3}/{4}\rho g\eta$ . The specific form of $P_b$ in the bulge-induced surface pressure, which is converted from the energy dissipation rate, is given as

(3.13) \begin{equation} \displaystyle {{P_b = \tfrac {3}{4}\rho g\eta }}. \end{equation}

As we know, the 2-D solitary wave solution governed by the KdV equation maintains a constant wave height, which is the result of the balance between nonlinear effects and dispersive effects. In the 2-D axisymmetric condition, however, the presence of the term ( ${{\tilde u}}/{{\tilde r}}$ ) in continuity equation disrupts the balance of dispersion and nonlinearity in the solitary wave, resulting in the decay of the leading wave and changing the wave shape. In the propagating process of bulge evolution, the bulge-induced term in (3.9) serves as an additional cause for the decay of the leading wave amplitude. Therefore, during the propagation process, the decay in the leading wave amplitude will exhibit two components. The first component is caused by the mass conservation, which signifies the decay resulting from the gradually increasing length in the round-direction. The second component is the bulge-induced term, representing the decay caused by energy transfer from the leading wave to the tail wave train.

3.2.2. Order-of-magnitude analysis

It should be noted that the effects of viscosity dissipation are neglected in the theoretical analysis. Unlike the diffusion process, the propagating process involves a relatively high-Re process, which means that viscous dissipation near the bottom has little influence on the wave shape (Miles Reference Miles1977). To show the effect of viscosity on the entire process, the ratio of the energy of the viscous term ( $E_v$ ) to the energy of bulge-induced pressure term ( $E_b$ ) is examined later. Additionally, according to Longuet-Higgins (Reference Longuet-Higgins1957), the boundary layer near the free surface may be considered negligible in the long wave approximation for high Re. Therefore, in the following analysis of wave height decay during the wave propagation process, the influence of viscosity dissipation, being a minor factor, is neglected.

Similar to the analysis of the diffusion process, to obtain a simple equation for the surface elevation, the assumption of weak nonlinearity and long wavelength is also required. However, when dealing with concentric waves, nonlinear analysis becomes more complex and some stretched variables are introduced. For radial variable, $r$ , a spreading parameter measured by $\lambda /r$ is introduced. Due to the radial invariance of their relative magnitudes in a channel flow with a slowly varying width, there exists a conserved quantity, $\eta /h \sim (r/h)^{-2/3}$ (Lamb Reference Lamb1932). Therefore, $\tilde {r}$ is stretched with $\varepsilon ^{3/2}$ . Moreover, time $t$ , is logically scaled by $r/c_0$ , where a characteristic time scale is naturally defined as the characteristic period of the linear waves. So, $\tilde {t}$ is also stretched with $\varepsilon ^{3/2}$ . Since the motion is driven by the free surface, which is very gentle, the characteristic fluid velocity $u$ , is also very gentle compared with the phase speed of the linear wave. Consequently, we choose $u$ to be of the order $O(\varepsilon )$ , resulting in the stretching of $\hat u = {\varepsilon ^{ - 1}}\tilde u$ . The fluid velocity $\tilde {w}$ is naturally stretched by $\varepsilon ^{-3/2}$ , as indicated by the continuity equation. Furthermore, the wave height $\tilde {\eta }$ is stretched by $\varepsilon ^{-1}$ , and the pressure $\tilde {p} = \tilde {\eta }$ is stretched by $\varepsilon ^{-1}$ naturally. The same stretching applies to the bulge-induced surface pressure. Therefore, stretched variables can be summarised as $\hat r = {\varepsilon ^{3/2}}\tilde r;\;\hat z = \tilde z;\;\hat t = {\varepsilon ^{3/2}}\tilde t;\;\hat u = {\varepsilon ^{ - 1}}\tilde u;$ $\hat w = {\varepsilon ^{ - 3/2}}\tilde w;\;\hat p = {\varepsilon ^{ - 1}}\tilde p;\;{\hat P_b} = {\varepsilon ^{ - 1}}{\tilde P_b};\hat \eta = {\varepsilon ^{ - 1}}\tilde \eta$ . For the long-time simulations of an outward propagating wave, it is desirable to switch to a moving reference frame by making the change of variables $(\hat \xi ,\hat \tau ) = ({\varepsilon ^{ - 1}}(\hat r - \hat t),\hat t)$ , so that the observer moves outwards in the radial direction with the linear long-wave speed. Similar to the diffusion process, appropriate scale transformations are necessary to ensure that the Re term is minimal compared with the convective term, making the effect of the viscosity term negligible. Additionally, ensuring that the order of the bulge-induced term is consistent with the convective term, allowing the impact of the nonlinear wave train on the leading wave to be clearly observed, yields $Re^* = {\varepsilon ^{3/2}}{{Re}};\ \hat P_b^* = {\varepsilon ^{1/2}}{{\hat P}_b}$ . With the help of $\varepsilon$ , all variables are appropriately stretched and of order $O(1)$ . Such stretched transformations can facilitate clear expansion of these stretched variables according to their order of magnitude $O(\varepsilon ^n)$ .

Expressed in terms of these newly introduced variables, the governing equations and boundary conditions can be written with only the leading order of $\varepsilon$ explicitly considered as

(3.14) \begin{equation} \left \{ \begin{array}{l} \displaystyle {\frac {\varepsilon }{{\hat \tau }}\hat u + \frac {{\partial \hat u}}{{\partial \hat \xi }} + \frac {{\partial \hat w}}{{\partial \hat z}} = 0}\\[13pt] \displaystyle {\varepsilon \frac {{\partial \hat u}}{{\partial \hat \tau }} - \frac {{\partial \hat u}}{{\partial \hat \xi }} + \varepsilon \hat u\frac {{\partial \hat u}}{{\partial \hat \xi }} + \varepsilon \hat w\frac {{\partial \hat u}}{{\partial \hat z}} = - \frac {{\partial \hat p}}{{\partial \hat \xi }} }\\[13pt] \quad \displaystyle {+ \frac {1}{{{\mathrm{R}}{{\mathrm{e}}^*}}}\left ( {\varepsilon \frac {{{\partial ^2}\hat u}}{{\partial {{\hat z}^2}}}} \right ) - {{\hat P}_b}^*\varepsilon \frac {1}{{\hat \eta }}{\left(\frac {{\partial \hat \eta }}{{\partial \hat \xi }}\right)^2} - {{\hat P}_b}^*\varepsilon \frac {{{\partial ^2}\hat \eta }}{{\partial {{\hat \xi }^2}}}} \\[13pt] \quad\displaystyle { - \varepsilon \frac {{\partial \hat w}}{{\partial \hat \xi }} = - \frac {{\partial \hat p}}{{\partial \hat z}}} - {\varepsilon ^{ - 1}} \end{array} \right .\quad \left \{ \begin{array}{l} \hat u = \hat w = 0\quad \mathrm{at}\quad \hat z = - 1 \\[5pt] \displaystyle {\hat w = \varepsilon \frac {{\partial \hat \eta }}{{\partial \hat \tau }} - \frac {{\partial \hat \eta }}{{\partial \hat \xi }} + \varepsilon \hat u\frac {{\partial \hat \eta }}{{\partial \hat \xi }}} \\[13pt] \quad\mathrm{at}\quad \hat z = \varepsilon \hat \eta (\hat \xi ,\hat \tau ) \\[6pt] \displaystyle {{\frac {{\partial \hat u}}{{\partial \hat z}}} = 0\quad \mathrm{at}\quad \hat z = \varepsilon \hat \eta (\hat \xi ,\hat \tau ) } \\[10pt] {\hat p_\infty } - \hat p = 0\quad \mathrm{at}\quad \hat z = \varepsilon \hat \eta (\hat \xi ,\hat \tau ). \end{array} \right . \end{equation}

3.2.3. Perturbation expansion

By introducing weakly nonlinear expansions (3.3) and free surface Taylor expansion (3.4) into (3.14), a straightforward equation for the surface elevation can be derived by reducing multiple-scale expansion equations to their $O(\varepsilon )$ terms as

(3.15) \begin{equation} \frac {{{{\hat \eta }_0}}}{{2\hat \tau }} + \frac {{\partial {{\hat \eta }_0}}}{{\partial \hat \tau }} + \frac {3}{2}{\hat \eta _0}\frac {{\partial {{\hat \eta }_0}}}{{\partial \hat \xi }} + \frac {1}{6}\frac {{{\partial ^3}{{\hat \eta }_0}}}{{\partial {{\hat \xi }^3}}} + \frac {3}{8}{\left(\frac {{\partial {{\hat \eta }_0}}}{{\partial \hat \xi }}\right)^2} + \frac {3}{8}{\hat \eta _0}\frac {{{\partial ^2}{{\hat \eta }_0}}}{{\partial {{\hat \xi }^2}}} = 0. \end{equation}

Currently, no known exact solutions exist for this type of variable-coefficient KdV equation. However, we can solve the equation by assuming that the solution adheres to the form of a solitary-wave solution of the KdV equation (Mandi, Roy & Chatterjee Reference Mandi, Roy and Chatterjee2020), then by using the conserved quantity, $\int _{ - \infty }^\infty {\hat \eta _0^2{\rm d}\xi }$ , the solution for the wave height and time can be derived as

(3.16) \begin{equation} {\hat \eta _0}^{ - 2} = \underbrace {(36/35)\hat \tau }_{\scriptstyle {\mathrm{bulge\hbox{-}induced\thinspace pressure }}\atop \scriptstyle {\mathrm{impact\thinspace term}}} + \underbrace {{C}{{\hat \tau }^{4/3}}}_{\scriptstyle {\mathrm{mass\thinspace conservation}}\atop \scriptstyle {\mathrm{impact\thinspace term}}}, \end{equation}

where $C$ is a constant. A detailed derivation is provided in Appendix B. It should be noted that this method is only applicable to slowly varying waves. In (3.16), the first term on the right-hand side originates from the contribution of the bulge-induced pressure, while the second term accounts for the contribution of the two-dimensional axisymmetric radial variation. When the influence of the bulge-induced pressure is neglected, (3.16) reduces to the classical C-KdV wave height decay scaling law.

3.2.4. Viscous effect estimation

Now, the influence of viscous dissipation during the wave propagation process will be analysed. For the sake of simplicity in analysis, the 2-D case is considered here. Ott & Sudan (Reference Ott and Sudan1969) presented the KdV equation considering the viscous dissipation given in (3.17) as

(3.17) \begin{equation} {\hat \eta _{0,\hat \tau }} + \frac {3}{2}{\hat \eta _0}{\hat \eta _{0,\hat \xi }} + \frac {1}{6}{\hat \eta _{0,\hat \xi \hat \xi \hat \xi }} = \frac {1}{{4\sqrt {\pi Re} }}\int \limits _{ - \infty }^\infty {\frac {{\partial {{\hat \eta }_0}}}{{\partial \hat \xi '}}\frac {{1 - {\mathrm{sgn}}\big( {\hat \xi - \hat \xi '} \big)}}{{\sqrt {\big | {\hat \xi - \hat \xi '} \big|} }}{\rm d}\hat \xi '}. \end{equation}

Similar to the derivation of the wave equation in the above 2-D axisymmetric case, the bulge-induced pressure in the 2-D propagation process takes the form $P_b = 0.25\rho g \eta ^3 / h^2$ . The wave equation in the 2-D case is shown in (3.18) as

(3.18) \begin{equation} \frac {{\partial {{\hat \eta }_0}}}{{\partial \hat \tau }} + \frac {3}{2}{{\hat \eta }_0}\frac {{\partial {{\hat \eta }_0}}}{{\partial \hat \xi }} + \frac {1}{6}\frac {{{\partial ^3}{{\hat \eta }_0}}}{{\partial {{\hat \xi }^3}}} = - \frac {3}{8}{{\hat \eta }_0}^2{\left(\frac {{\partial {{\hat \eta }_0}}}{{\partial \hat \xi }}\right)^2} - \frac {1}{8}{{\hat \eta }_0}^3\frac {{{\partial ^2}{{\hat \eta }_0}}}{{\partial {{\hat \xi }^2}}}. \end{equation}

To assess the impact of viscous effects on the propagation process, we independently integrate the terms involving wave height ( ${{\hat \eta }_0}$ ) multiplied by the right-hand sides of (3.17) and (3.18). Due to the complexity of the right-hand side of (3.17), the equation can be solved by using the Fourier transform. Consequently, we derive the energy contributions from viscosity and bulge-induced terms, as demonstrated in (3.19), where ${\mathord {\buildrel {\lower 3pt\hbox{$\scriptscriptstyle \frown $}} \over \eta } _0}(k,\hat \tau )$ is the Fourier transform of ${{{\hat \eta }_0}}(\hat \xi ,\hat \tau )$ . Here, ${\mathord {\buildrel {\lower 3pt\hbox{$\scriptscriptstyle \frown $}} \over \eta } _0}(-k,\hat \tau )$ is the complex conjugate to the ${\mathord {\buildrel {\lower 3pt\hbox{$\scriptscriptstyle \frown $}} \over \eta } _0}(k,\hat \tau )$ with the wavenumber $k$ . By substituting the expression of solitary wave solution based on Boussinesq theory in the spectral domain and the space domain into (3.19), doing the integration, it can be derived that ${E_v}/{E_p} \sim 1/\sqrt {{Re}} \ll 1$ . Hence, neglecting the viscous dissipation is reasonable in the propagation process.

(3.19) \begin{equation} {E_v} = - \frac {{\pi 2\sqrt 2 }}{{\sqrt {Re} }}\int \limits _0^\infty {\sqrt {\left | k \right |} {{\mathord {\buildrel {\lower 3pt\hbox{$\scriptscriptstyle \frown $}} \over \eta } }_0}\left ( { - k,\hat \tau } \right ){{\mathord {\buildrel {\lower 3pt\hbox{$\scriptscriptstyle \frown $}} \over \eta } }_0}\left ( {k,\hat \tau } \right ){\rm d}k} ;\quad {E_p} = - \int \limits _{ - \infty }^\infty {\frac {1}{8}{{\hat \eta }_0}^3{{\left(\frac {{\partial {{\hat \eta }_0}}}{{\partial \hat \xi }}\right)}^2}} d\hat \xi .\end{equation}

Figure 9.

Verification of scaling laws between wave height and time of propagation process. All data points have been translated to start from an initial value of zero. The subscript $i$ denotes the selected initial value. The numbers following $h$ , $H$ , $W$ represent their values, expressed in SI units. The black solid lines, each with a slope of 45 degrees, indicate the trend given by the theoretically derived scaling law (‘soli’, ‘cos’, and ‘cyl’ represent the initial bulge profile as a solitary, cosine wave and cylindrical shapes, respectively).

3.2.5. Scaling law discussion

Figure 9 presents the variations in wave height under different conditions obtained through DNS. The conditions of multiple stationary water depths and initial bulge heights are set. In addition to using the solitary wave profile for the initial bulge, other initial profiles such as cosine wave and cylindrical shape are also tried. The previously derived scaling law for wave height is examined. Here, $C$ is a constant, its value calculated as the average of all $C$ values obtained by substituting all data points into (3.16) in each condition. It can be observed that the numerical results align well with the wave height decay over time relationship in (3.16), demonstrating that the modelling and analysis of the wave propagation process are reasonable. Moreover, to capture the limitations of this theoretical model and identify the applicability boundaries of the derived scaling law, the numerical simulations in this study further examined conditions that deviate significantly from the original assumptions, specifically those involving stronger nonlinearities with $ H / h \gt 0.1$ . To further refine the specific applicability range, we have introduced expanded extreme-case simulations to determine the limitations of the scaling law, as shown in figure 9. We find that the derived scaling law is not strictly confined to weakly nonlinear conditions, but also remains valid in cases with moderately strong nonlinearity. However, when the initial wave height becomes greater than the water depth $H/h\gt 1$ , deviations from the theoretical scaling law begin to emerge. Therefore, the wave height attenuation scaling law derived in this study is applicable to shallow water wave propagation, when the initial wave height is smaller than the water depth with a moderately strong nonlinearity.

It is worth noting that the $C$ values of the mass conservation impact term shown in figure 9 range approximately between $O(10^0)$ and $O(10^1)$ , which are relatively larger than the coefficient $36/35$ of the bulge-induced pressure impact term. This suggests that, at smaller time scales, the bulge-induced pressure impact term has some influence on the decay of the leading wave. However, on larger time scales, the decay of the leading wave is primarily governed by the mass conservation impact term. In the C-KdV equation framework, the attenuation of the leading wave is fundamentally governed by radial expansion necessary to maintain mass conservation. This physical mechanism establishes the coefficient $C$ in the wave height scaling law as a critical indicator of radial mass transport efficiency within the leading wave. Specifically, enhanced $C$ values correspond to accelerated mass redistribution along the radial direction, thereby intensifying wave height decay rates. The magnitude of $C$ arises from two competing factors: the spatial extent of the leading wave and its phase velocity. Increasing wave dimensions suppress $C$ values through enhancing the cross-sectional area, while greater phase velocities amplify $C$ by facilitating rapid mass redistribution. Additionally, considering the zero-order phase velocity $c_0 = \sqrt {gh}$ and substituting $t = r/c_0$ into (3.16), the decay relationship between ${\hat \eta }_0$ and $\hat {r}$ can be obtained. This relationship has the same form as (3.16), but with different coefficients.