The lifecycle of axisymmetric internal solitary waves

The generation and evolution of solitary waves by intrusive gravity currents in an approximate two-layer fluid with equal upperand lower-layer depths is examined in a cylindrical geometry by way of theory and numerical simulations. The study is limited to vertically symmetric cases in which the density of the intruding fluid is equal to the average density of the ambient. We show that even though the head height of the intrusion decreases, it propagates at a constant speed well beyond 3 lock radii. This is because the strong stratification at the interface supports the formation of a mode-2 solitary wave that surrounds the intrusion head and carries it outwards at a constant speed. The wave and intrusion propagate faster than a linear long wave; therefore, there is strong supporting evidence that the wave is indeed nonlinear. Rectilinear Korteweg-de Vries theory is extended to allow the wave amplitude to decay as r with p= 2 and the theory is compared to the observed waves to demonstrate that the width of the wave scales with its amplitude. After propagating beyond 7 lock radii the intrusion runs out of fluid. Thereafter, the wave continues to spread radially at a constant speed, however, the amplitude decreases sufficiently so that linear dispersion dominates and the amplitude decays with distance asr−1.


Introduction
A gravity current arises when a fluid of one density flows horizontally into an ambient fluid of a different, or varying, density.In a stratified fluid, a gravity current can propagate at an intermediate level between the upper and lower boundaries of the ambient along its level of neutral buoyancy.For this reason, these currents are specifically referred to as "intrusive gravity currents" or "intrusions".
Correspondence to: B. R. Sutherland (bsuther@ualberta.ca)Laboratory experiments have been performed to examine the simplest circumstance of "symmetric" intrusions that propagate in an approximate two-layer fluid with equal upper-and lower-layer depths in which the intrusion density was the average ambient density (Britter and Simpson, 1981;Faust and Plate, 1984;Rooij et al., 1999;Lowe et al., 2002).By symmetry, such studies are equivalent to the examination of surface gravity currents moving across a free-slip boundary with a stratification of the ambient near the surface, an idealization of the dynamics governing river plumes (Nash et al., 2009).Unlike gravity currents in uniform ambients, which are predicted to decelerate after propagating 6 to 10 lock lengths (Rottman and Simpson, 1983), symmetric intrusions at sufficiently thick interfaces were found to propagate well beyond this distance at a constant speed (Mehta et al., 2002;Sutherland and Nault, 2007).In a rectilinear geometry, intrusions maintained a constant speed beyond 22 lock lengths.For small, but finite interface thicknesses, the measured speeds were faster than the linear long wave speed, suggesting that the intrusions first created a closed-core solitary wave and then were carried by the wave.
The dynamics of radially spreading gravity currents are qualitatively different because energy and mass conservation require the head height to decrease and hence the current should decelerate soon after release due to the reduction of the horizontal pressure gradient driving the flow.Lockrelease experiments (Huppert and Simpson, 1980;Didden and Maxworthy, 1982;Huq, 1996;Hallworth et al., 1996;Patterson et al., 2006) demonstrated that axisymmetric bottom-propagating gravity currents maintain a constant speed up to 3 lock radii.Thereafter, in the self-similar regime, the position of the front changed as t 1/2 , consistent with the prediction of shallow water wave theory (Huppert and Simpson, 1980).By symmetry, Zemach and Ungarish (2007) used shallow water theory to extend this prediction to axisymmetric intrusion propagation, again predicting that the intrusion decelerates shortly after release from the lock.
Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
Contrary to this prediction, however, full-depth lock release experiments demonstrated that axisymmetric intrusions propagate at a constant speed well beyond 3 lock radii (Sutherland and Nault, 2007), suggesting that, as in a rectilinear geometry, solitary waves were responsible for transporting the lock-fluid long distances at a constant speed.This observation puts into question the applicability of shallow water theory, which necessarily filters nonhydrostatic dynamics.
The theory for internal solitary waves in a cylindrical geometry is limited.In a rectilinear geometry, the weakly nonlinear Korteweg-de Vries (KdV) model is appropriate if nonlinearity and nonhydrostatic dispersion are comparable and small.Adapting the fully nonlinear Dubreil-Jacotin-Long (DJL) equation (Dubriel-Jacotin, 1937;Long, 1953Long, , 1956)), the large-amplitude structure of steady rectilinear closed-core (Davis and Acrivos, 1967;Tung et al., 1982;Brown and Christie, 1998) and leaky (Derzho and Grimshaw, 2007) solitary waves have been modelled.In particular, White and Helfrich (2008) used DJL theory to examine the propagation of gravity currents and intrusions in a rectilinear geometry.Because this theory assumes the system is in steady state, it is unclear that it can be straightforwardly adapted to a cylindrical geometry in which the wave amplitude must decrease with radius and hence in time.
Following up upon the work of Miles (1978), Weidman and Velarde (1992) developed a theory for axisymmetric internal solitary waves within a stratified ambient fluid.They predicted an amplitude dependence upon distance as r −2/3 .But, in an earlier paper Weidman and Zakhem (1988) compared this amplitude prediction with experiments by Maxworthy (1980) and found that the amplitude differed by up to 34%.They explained that their weakly nonlinear theory did not capture the experimental results because the trapping of the mixed fluid was "a clear manifestation of strong nonlinear effects".They admitted, however, that their equation was asymptotically inconsistent and so its predictions were not necessarily reliable.The theory presented by Weidman and Velarde (1992) likewise relies upon inconsistent assumptions.
In this work, we perform the first detailed numerical simulations of radially spreading axisymmetric intrusions and solitary waves, the results of which are tested against laboratory experiments and a heuristic theory that is a straightforward adaptation of rectilinear KdV theory.The fact that the behaviour observed in our simulations more closely resembles the heuristic prediction and that it differs from the behaviour predicted by the cylindrical KdV equation puts into question the applicability of the latter to solitary waves generated by radial intrusions.
The paper is organized as follows.In Sect.2, the theoretical speed of a vertically symmetric intrusion is given as a function of interface thickness and the theories of axisymmetric linear waves and rectilinear solitary waves are also reviewed.Here we present a heuristic theory for cylindrical solitary waves and compare their corresponding equation with the KdV cylindrical solitary wave equation whose derivation is given in detail in the Appendix.The formulation of the fully nonlinear simulations and their results are presented in Sect.3. Conclusions regarding the comparison of simulations to weakly nonlinear theories are given in Sect. 4.

Theory
This work involves the study of radially spreading intrusions and the interfacial waves they generate.We begin by adapting existing theory for rectilinear gravity currents in uniform density fluid and in uniformly stratified fluid to develop a prediction for the speed of intrusions at interfaces of arbitrary thickness.This analysis is restricted to the study of symmetric intrusions, meaning that the density of the intrusion is the average ambient density and the background density gradient is itself symmetric in the vertical.In Sect.3.2, we proceed to compare the observed speed of radial intrusions with the rectilinear theory prediction and so evaluate the effect of initial curvature upon setting the radial intrusion speed.Being symmetric, the intrusion can most efficiently excite a mode-2 interfacial wave, which bulges upwards above and downwards below the mid-depth of the interface.We briefly review the theory for small-amplitude long interfacial waves with this mode-2 structure.Finally, we review Korteweg-de Vries (KdV) theory for rectilinear solitary waves in continuously stratified fluid and then describe an adaptation of this theory for radially spreading solitary waves, which is later compared to observations.

Rectilinear speed of a symmetric intrusion at a thick interface
The speed of a rectilinear intrusion at an interface of finite thickness was first predicted by Ungarish (2005) using a shallow water model, although that theory underpredicted the experimental speeds measured by Faust and Plate (1984).White and Helfrich (2008) adapted DJL theory to arrive at a fully nonlinear description of intrusions in stratified ambients.Although successful for moderately thick interfaces, the theory unphysically required the head height of the intrusion in a linearly stratified ambient to tend to zero.In this limit, the theory overpredicted the intrusion speed measured in experiments and simulations by Bolster et al. (2008).
Here, we take the simplest approach to predict the intrusion speed while ensuring that our result is in agreement with the predictions for two-layer and linearly stratified ambients.The background density profile, ρ, is assumed to have the piecewise-linear form, where h is the thickness of the interface and ρ i = 1 2 (ρ U + ρ L ) is the average of the upper and lower layer densities.To take advantage of symmetry, the level z = 0 is positioned at the mid-depth of the fluid, which has a total depth H .The intrusion can be viewed as a symmetric expansion of a surface-or bottom-propagating gravity current within an ambient of total depth h T = H 2 , as illustrated in Fig. 1.We first formulate the speed of a gravity current with density ρ i moving in an ambient with density given by Eq. ( 1) Benjamin (1968) predicted that an energy conserving gravity current released from a lock of height h T will propagate into a uniform ambient with a head height of h T /2.Assuming the head height changes negligibly in a stratified ambient, the mean density of the ambient over the depth of the gravity current is given by where δ h = h H .The gravity current speed, U gc , given by Benjamin (1968) and Ungarish (2006) is By extension, the speed of an intrusion, U i , moving along the interface of a stratified fluid with a density profile given by Eq. ( 1) is predicted by combining Eqs. ( 2) and (3) and letting h T = H /2 to give where U 0 is the speed of a symmetric intrusion in a two-layer fluid The reduced gravity, g LU , is given by where ρ 00 is a characteristic density.In deriving, Eq. ( 4) we have assumed that the intrusion speed is set by the density difference between the mean density of the ambient over the head height of H /2. As illustrated by the dashed line in Fig. 2, this intrusion speed, U i , is predicted to decrease as the relative thickness of the interface, δ h , increases.This formula correctly predicts the rectilinear intrusion speed in the two-layer limit (δ h → 0) and in the uniformly stratified limit (δ h → 1).2) are compared to the predicted rectilinear speed (dashed line) defined by Eq. ( 4) and the linear long wave speed (solid line) defined by Eq. ( 11).Experimental speeds, measured by Sutherland and Nault (2007), are plotted as circles.The largest experimental error bars are also shown in the lower right hand corner.

Long small-amplitude waves at a thick interface
Symmetric intrusions may efficiently excite solitary waves if they travel faster than long mode-2 interfacial waves at a thick interface.Rectilinear intrusion speeds were thus compared with long mode-2 waves in the x-z plane with a piecewise-constant three-layer density profile (Mehta et al., 2002;Sutherland and Nault, 2007).Here we find the speed of long axisymmetric interfacial waves in a threelayer fluid with piecewise-linear background density given by Eq. ( 1).By applying no-normal-flow boundary conditions at z = ± H 2 , bounded solutions of the streamfunction have the form ψ(r,z,t) = A ψ J 1 (kr)φ(z)e −iωt , where A ψ is the amplitude, k is the radial wavenumber, ω is the wave frequency, J 1 is the Here the vertical wavenumber of disturbances within the interface is in which the squared buoyancy frequency of the thick interface is Using ( 8), k and ω are implicitly related by the dispersion relation In the long wave limit, this becomes tan where c 0 = ω/k is the long wave speed.This speed is plotted by the solid line in Fig. 2. The same result would be found for rectilinear interfacial waves: the structure of the waves changes, but the long wave speed is the same.Comparing this result to Eq. ( 4), we find that the intrusion speed is supercritical (U i > c 0 ) if δ h 0.5 and subcritical otherwise.Hence, if radially spreading intrusions propagate at speeds comparable to those of rectilinear intrusions, we expect solitary waves will be excited if the intrusion moves along a sufficiently thin interface (δ h 0.5) but not so thin that a mode-2 wave cannot be established by the widening of the interface (δ h > 0).It should be noted that fully nonlinear waves can exist in a band around c 0 , as described by White and Helfrich (2008).

Internal solitary waves
In a rectilinear geometry, the KdV equation is widely used to model internal solitary waves in a stratified fluid.The formulation by Benney (1966) predicts that, in the Boussinesq approximation, the vertical displacement field, ξ , has the form in which a(x,t) = a(x − ct) satisfies the KdV equation that includes the advection term: Here the subscripts denote derivatives and the constants γ and β are given by (Benney, 1966) where c 0 is the linear long wave speed.For a stratified ambient with the density profile given by Eq. ( 1), we exploit symmetry and consider only the upper half of the domain (h T = H 2 ).The corresponding vertical structure function, φ, is given by the long wave limit of Eq. ( 7) with 0 ≤ z ≤ H 2 and the long wave speed, c 0 , satisfies Eq. ( 11).
The solution of Eq. ( 13) which assumes an isolated, steadily propagating disturbance is The speed, c, and width, λ, of the wave depend on the maximum displacement amplitude, a 0 , by The KdV coefficients along with the wave speed and width are shown in Fig. 3 as functions of interface thickness.
From the numerical simulations, it is easier to diagnose the structure of the vertical velocity rather than the vertical displacement field.The vertical velocity derived from the vertical displacement field is in which the coefficient satisfies where A w is defined as the maximum value of w.In these expressions, the value 2.598 is determined from the empirically computed maximum value of sech 2 tanh.The challenge is to extend these results to describe a solitary wave that spreads radially.For a solitary wave spreading as an expanding ring with sufficiently large radius, we may assume that the front on any point along the circumference negligibly feels the curvature of the front but that the amplitude nonetheless decreases with radius r as r −1/2 , a result of energy conservation.Thus, by extension of Eq. ( 17), we predict that the vertical velocity field should evolve as in which the amplitude A w0 is measured at a distance r 0 where the axisymmetric solitary wave is first generated.
For measured A w0 , we can use Eqs.( 16) and ( 18) to determine the initial displacement amplitude, a 0 , the solitary wave speed, c, and width, λ.As the wave spreads radially, A w = A w0 (r 0 /r) 1/2 decreases and these equations predict that the displacement amplitude and speed should decrease and its width should increase.Eventually, the amplitude should be so small that weakly nonlinear effects associated with wave steepening should become so small that linear dispersion dominates the evolution.
In a more rigorous approach, we attempted to adapt the theory of Weidman and Velarde (1992) for axisymmetric internal solitary waves to include the wave amplitude decay as r −1/2 .The method, as outlined in Appendix A, results in the same leading-order equation as Weidman and Velarde (1992): in which γ and β are given by Eq. ( 14).Consistent with the observation by Weidman and Zakhem (1988) concerning the analogous result by Miles (1978) for surface waves, the cylindrical solitary wave equation is asymptotically inconsistent in that Eq. ( 20) is valid only for r r 0 in which case the balance between nonlinear steepening and dispersion may no longer be valid.
By contrast, the heuristic solution of the form Eq. ( 19) is established for r r 0 .The corresponding vertical displacement field of the form a(r,t) = a 0 (r 0 /r) 1/2 sech 2 [(r − ct)/λ] does not satisfy Eq. ( 20) but rather, through manipulation of rectilinear KdV theory, satisfies the leading-order equation Compared to Eq. ( 20), here the nonlinear term is of order (r/r 0 ) 1/2 , which seems to suggest that the balance between nonlinearity and dispersion exists for large r.However, because a ∼ r −1/2 , the nonlinear term will eventually become negligible and the motion will then be governed by This equation simply describes a dispersing linear wave and using the method of stationary phase it can be shown that in this regime, the equally competing effects of geometrical spreading and dispersion cause the amplitude of a long wave to decrease as r −1 .
Through comparison of these theories with the following results of numerical simulations, we will show that the heuristic solution, Eq. ( 19), more accurately represents the life-cycle of axisymmetric solitary waves generated by intrusions.Thus we show that the admittedly inconsistently derived cylindrical solitary wave equation of Weidman and Velarde (1992) and our extension of it in Appendix A, indeed provide an inaccurate description of axisymmetric solitary waves.However it must be emphasized that Eq. ( 21) is heuristic: whereas the balance in rectilinear KdV theory for a disturbance of length L is aL 2 ∼ 1, the nonlinear dispersive "balance" is aL 2 ∼ r −1/2 , and so the nature of competition between nonlinearity and dispersion changes with radial distance.

Description of code
Fully nonlinear numerical simulations were performed to examine the evolution of axisymmetric solitary waves in a two-layer Boussinesq fluid with finite interface thickness.
The code solves the cylindrical Navier-Stokes equations with the assumption that the azimuthal velocity is zero everywhere.The resulting two coupled partial differential equations for the azimuthal vorticity, ζ , and perturbation density, ρ, fields are given by where ∇ 2 is the Laplacian in cylindrical coordinates.The radial and vertical components of the velocity field are given by u and w, respectively.Due to the axisymmetric geometry, a streamfunction, ψ, can be defined such that u = − ∂ψ ∂z and w = 1 r ∂ψ ∂r .This is related implicitly to the vorticity field by Given ζ at a particular time, Fourier-Bessel transforms are used to invert Eq. ( 25) to find ψ.From this, u and w are computed and then Eqs. ( 23) and ( 24) are used to advect the vorticity and density fields.At each time step, a passive tracer field is also advected.
At the free-slip boundaries, the no-normal flow condition, u • n = 0, is imposed and it is assumed that ζ = 0.The numerical code approximates spatial derivatives using secondorder finite difference methods.The evolution equations are then stepped forward in time using a leap-frog scheme and, to minimize time splitting errors, an Euler backstep is taken every 20 time-steps.
To model the full-depth lock release experiments conducted by Sutherland and Nault (2007), a domain with a radius of R = 45 cm and a height of H = 10 cm was used.The equations were solved on a staggered grid with 1025 radial points and 257 vertical points.To examine the long-time evolution of the system, simulations with a domain radius of R = 80 cm were also performed.
The code was initialized by a density field that mimicked the initial density of the experimental setup, given by where r 0 =6 cm is the radius of the lock.Simulations were completed with both a piecewise linear background density, ρ(z), given by Eq. ( 1) and a smooth hyperbolic tangent profile given by, For all simulations, the density within the lock was the average ambient density, ρ i = 1 2 (ρ L + ρ U ).The concentration of the passive tracer field was initially set to unity over 0 < r < r 0 and zero elsewhere.
The physical parameters used in the simulations were g = 980.6 cm/s 2 and ρ 00 =1.0 g/cm 3 .To prevent the code from becoming numerically unstable while maintaining a reasonable computation speed, a spatially varying piecewiselinear viscosity was prescribed in the simulations.To damp out small-scale noise created by the collapse of the lock fluid, a viscosity of 0.1 cm 2 /s was used where 0 < r < 2r 0 .Further away from the lock, where r > 3r 0 , the physical value of ν = 0.01 cm 2 /s was used.Between these regions the viscosity varied linearly with r.The dynamics of the flow were unaffected using this viscosity because the Reynolds numbers were on the order of 10 3 which suggests that viscosity did not govern the dominant motion.To confirm this, a high resolution simulation with a uniform viscosity of 0.01 cm 2 /s was performed, this taking many days rather than many hours to run.The propagation and speed of the intrusion near the lock remained unchanged.The diffusivity of salt water, κ, was set to be everywhere equal to ν, even though its physical value is 10 −5 cm 2 /s.The purpose of doing this is again for numerical stability.Nonetheless the value of κ was sufficiently small that molecular diffusion had a negligible influence on the flow.

Results
In agreement with the experimental observations by Sutherland and Nault (2007), axisymmetric intrusions at a thin interface (δ h ≤ 0.2) were observed to propagate beyond 8r 0 at a constant speed as shown in Fig. 4. On the contrary, axisymmetric gravity currents have been observed to decelerate as early as 4r 0 (Huppert and Simpson, 1980;Patterson et al., 2006), therefore, it is believed that the stratification of the thin interface allows for the formation of a wave.The simulations show that intially the wave and intrusion propagate outwards together, however, when the intrusion decelerates to a stop, the wave continues to spread radially.As shown in Fig. 4, this deceleration, and hence separation, was observed to occur at smaller r for increasing δ h .The location of the intrusion front versus time for simulations with ρ L =1.0530 g/cm 3 and ρ(z) given by Eq. ( 1).At a thick interface (δ h ≥ 0.60), the intrusion begins to decelerate around 4r 0 , whereas at a thin interface (δ h ≤ 0.2) the intrusion maintains a constant speed beyond 8r 0 .
For a range of simulations, the initial intrusion speeds, C, were measured by calculating dr i /dt for 2r 0 ≤ r i ≤ 3r 0 .These speeds are plotted in Fig. 2 for intrusions in ambients with both piecewise linear and hyperbolic tangent profiles and the results are in agreement with those measured in laboratory experiments (solid circles).Consistent with the rectilinear theory outlined in Sect.2.1, as the thickness of the interface increases, the intrusion speed decreases.However, for all δ h , axisymmetric intrusions travel more slowly than the predicted speed in a rectilinear geometry, given by Eq. ( 4).This observation is consistent with axisymmetric gravity current experiments (Huppert and Simpson, 1980;Patterson et al., 2006) and simulations (Zhang et al., 2010), in which the observed front speeds were about 0.8U gc .Figure 2 also shows that compared to the long-wave speed, symmetric intrusions travel more quickly if δ h 0.4.
To examine the long-time evolution of the system, sideview snapshots are shown in Fig. 5 for a simulation with an interface thickness of δ h = 0.2 and an ambient density given by Eq. ( 27).The intrusion excites a mode-2 varicose wave which surrounds the intrusion head and carries it outwards at a constant speed.During the propagation, some lock fluid escapes rearward until the intrusion runs out of fluid and decelerates to a stop around 8r 0 .Beyond this distance, the wave continues to spread radially at a constant speed.This is illustrated by the dashed line in Fig. 6, where the wave's position, r w , is the location of maximum isopycnal displacement at a height of z/H = 0.15.
The evolution of the wave amplitude, A w , was determined by taking a radial time series of the vertical velocity field at a height of z/H = 0.15.To ignore interactions of the wave with the intrusion, the amplitude of the leading wave crest was measured.The results of this analysis are illustrated in Fig. 6 where it evident that the wave amplitude initially decays as r −1/2 w .This behaviour is anticipated on the basis of energy conservation for a non-dispersive wave.Beyond 8r 0 , the wave separates from the intrusion and its amplitude becomes sufficiently small that linear dispersion increases the decay rate to r −1 w .This late-time behaviour is consistent with the additional effects of linear dispersion on a small amplitude wave as discussed in Sect.2.3.
The observation that the wave initially travels faster than the long wave speed provides evidence that it is nonlinear upon generation.Because we found the amplitude does not decay as r −2/3 as predicted by Weidman and Velarde (1992) for a cylindrical solitary wave, here we take the simplest approach, of adapting the KdV model to include the observed amplitude decrease with radius as r −1/2 .
The results of the simulation shown in Fig. 5 were compared to the theory in Sect.2.3 by first estimating the height, h/2, of the theoretical interface.Because the simulation had a continuously stratified ambient, ρt , given by Eq. ( 27), the height at which the vertical velocity was a maximum (z/H 0.15) was assumed to correspond to z = h/2, where the normalized streamfunction amplitude satisfies φ(z) = 1.Radial slices of the vertical velocity field were then taken at several times, as shown in Fig. 7a.For each slice, the maximum amplitude, A w , was determined and Eqs. ( 16) and ( 18) were solved for a 0 , c and λ.
As shown in Fig. 7b, the profiles collapse onto a theoretical curve when the amplitude is scaled by r −1/2 w and the radial extent is shifted by 2r 0 + c t and scaled by λ.It should be noted that the reference location of 2r 0 was chosen to ignore the initial generation of the wave caused by the intrusion.For the intermediate times, t/t 0 = 5.0 and 7.5, the simulation results are in excellent agreement with the theory.However, at t/t 0 = 10.0, the wave amplitude has become small and linear dispersion has slightly increased the broadening of the wave.The slight discrepancy at t/t 0 = 2.5 can be explained by the strong interaction between the wave and intrusion at early times.

Conclusions
Laboratory experiments and numerical simulations show that axisymmetric intrusions propagate at a constant speed over distances much longer than 3 lock radii.Their speed is dependent upon the non-dimensional thickness of the interface, decreasing from 0.8U 0 to 0.5U 0 as δ h increases from 0 → 1.These speeds are up to 20% slower than those of intrusions in a rectilinear geometry.
For interfaces satisfying 0 < δ h 0.4, the intrusions were observed to excite a mode-2 varicose wave that surrounded the intrusion head.The wave then carried the intrusion along at a speed greater than the predicted linear long wave speed.Profiles of the wave just above the intrusion   5. Snapshots of the normalized vertical velocity field, w/w max , obtained from a numerical simulation of an axisymmetric intrusion in an ambient fluid with a background density, ρt , given by Eq. 27 with a non-dimensional interface thickness of δ h = 0.2.The thick black lines outline the intrusion profile at each time and illustrate that the intrusion head is being carried outward by a wave.In each plot, the vertical velocity field is normalized by the maximum amplitude of the wave at t/t 0 = 0.5, where t 0 = r 0 /U 0 .The colourbar limits are scaled by t/t 0 − 0.5, making it evident that the amplitude of the wave decays from its maximum value as t −1/2 .The wave is observed to propagate at a constant speed (i.e.r ∼ t); therefore, the amplitude of the wave is decreasing as r −1/2 as is predicted by linear theory.A reference amplitude, A w 0 , is defined such that A w = A w 0 (r w /r 0 ) −0.48 for 2r 0 ≤ r w ≤ 3r 0 .The circles correspond to t/t 0 = 2.5, 5.0, 7.5 and 10.0, for which wave profiles are illustrated in Fig. 7.The radial position of the wave versus time is indicated by the dashed line.A reference time, t * , is defined such r w = dr w dt (t − t * ) for 4r 0 ≤ r w ≤ 10r 0 .The curves are compared with lines of constant slope as indicated.
head were compared to a rectilinear solitary wave KdV theory adapted to a cylindrical geometry in which the wave amplitude decreased as r −1/2 .The wave profiles collapsed onto a theoretical curve suggesting that the amplitude decay was slow enough that it did not significantly affect the propagation of the wave.From a single measurement of wave amplitude and the assumption that the amplitude decays as r −1/2 , KdV theory was able to predict the amplitude, speed and spread of the wave during its nonlinear evolution phase after generation.
After propagating beyond 8r 0 , the intrusion ran out of fluid and the rate of amplitude decay of the wave increased to r −1 .Through a combination of weakening nonlinearity and increasing dominance of dispersion, the wave continued to propagate at a constant speed.
The fact that the solitary wave is observed to evolve closer to the heuristic prediction Eq. ( 19) than to the prediction of the cylindrical solitary wave equation suggests that a different mathematical approach should be taken in developing an appropriate weakly nonlinear theory.Finally, we note that the KdV solitary wave theory we have adapted does not account for the leaky closed-core behaviour associated with a wave carrying fluid associated with an intrusion.Corresponding normalized and shifted profiles, where t is the time elapsed since r w = 2r 0 , p = − 1 2 and A w 0 is a reference amplitude defined in Fig. 6.
The DJL equation, though used to examine intrusions and large-amplitude solitary waves in a rectilinear geometry (White and Helfrich, 2008), is less well adapted to study axisymmetric intrusions and waves because their dynamics are intrinsically unsteady and identical upstream and downstream boundary conditions cannot straightforwardly be assumed.In future work we will adapt this idealized study to more realistic oceanographic circumstances in order to predict the behaviour of radially spreading solitary waves generated by river plumes.

Fig. 1 .Fig. 2 .
Fig. 1.The domain and the density profile used to derive a theory for the speed of intrusions as a function of interface thickness.
Fig. 3. (a) Coefficients of the KdV equation as functions of interface thickness.The nonlinear coefficient, γ , is shown by the solid line and the dispersion coefficient, β, is shown by the dashed line.(b) The speed, c, (solid line) and width, λ, (dashed line) of a rectilinear solitary wave as a function of interface thickness.
Fig.4.The location of the intrusion front versus time for simulations with ρ L =1.0530 g/cm 3 and ρ(z) given by Eq. (1).At a thick interface (δ h ≥ 0.60), the intrusion begins to decelerate around 4r 0 , whereas at a thin interface (δ h ≤ 0.2) the intrusion maintains a constant speed beyond 8r 0 .
Fig.5.Snapshots of the normalized vertical velocity field, w/w max , obtained from a numerical simulation of an axisymmetric intrusion in an ambient fluid with a background density, ρt , given by Eq. 27 with a non-dimensional interface thickness of δ h = 0.2.The thick black lines outline the intrusion profile at each time and illustrate that the intrusion head is being carried outward by a wave.In each plot, the vertical velocity field is normalized by the maximum amplitude of the wave at t/t 0 = 0.5, where t 0 = r 0 /U 0 .The colourbar limits are scaled by t/t 0 − 0.5, making it evident that the amplitude of the wave decays from its maximum value as t −1/2 .The wave is observed to propagate at a constant speed (i.e.r ∼ t); therefore, the amplitude of the wave is decreasing as r −1/2 as is predicted by linear theory.

Fig. 6 .
Fig.6.The amplitude of vertical velocity, A w , versus radial position, r w , for the wave at a height of z/H = 0.15 (solid line).A reference amplitude, A w 0 , is defined such that A w = A w 0 (r w /r 0 ) −0.48 for 2r 0 ≤ r w ≤ 3r 0 .The circles correspond to t/t 0 = 2.5, 5.0, 7.5 and 10.0, for which wave profiles are illustrated in Fig.7.The radial position of the wave versus time is indicated by the dashed line.A reference time, t * , is defined such r w = Fig. 7. (a) Horizontal profiles of the vertical velocity field at a height z/H = 0.15.(b) Corresponding normalized and shifted profiles, wheret is the time elapsed since r w = 2r 0 , p = − 1 2 and A w 0 is a reference amplitude defined in Fig.6.