3.1. Steady drag-free flow

We seek fully time-dependent solutions to (3.1), (3.2) in which the drag is negligible ( $C_{D}=0$ ), but, as previously noted, will find that these time-dependent solutions are in non-similarity form, and much of the time-dependent solution is in fact steady. We therefore first look for steady solutions. The governing equations (3.1) and (3.2) can be written in steady-state form as

where (3.6b ) follows directly from a simplification of (3.2). Integration of these ordinary differential equations, which correspond to mass conservation and a Bernoulli equation, requires two boundary conditions. The constant associated with (3.6a ) is supplied by the condition (3.4) at source. The constant associated with (3.6b ) cannot be supplied by the front condition (3.5), since the inherently unsteady front condition cannot apply to a steady-state intrusion. We instead specify the constant of integration of (3.6b ), writing

where $E$ is for now an arbitrary positive constant, which we shall show is determined at the source of the intrusion.

With these constants determined, the solutions to (3.6a,b ) can be obtained algebraically (Baines Reference Baines2013). There are two solutions for the velocity and thickness fields: a supercritical solution in which $\mathit{Fr}=2|u_{r}|/h\geqslant 1$ ,

and a subcritical solution in which $\mathit{Fr}\leqslant 1$ ,

where

Both of these solutions exist only for $r\geqslant 1/E$ , and at $r=1/E$ the flow is critical, i.e. $\mathit{Fr}=1$ . We will show that it is always the supercritical solution (3.8) that is realised as part of the fully time-dependent motion, and plot this supercritical solution in figure 1, scaling $h$ and $u_{r}$ by $E^{-1/2}$ and $r$ by $E$ to remove the dependence on the parameter $E$ . In the supercritical solutions, $u_{r}$ tends to a constant and $h$ decays with increasing radial distance $r$ , meaning that the Froude number $\mathit{Fr}=2|u_{r}|/h$ diverges with increasing $r$ .

Figure 1.

The scaled velocity $u_{r}/\sqrt{E}$ (dashed line) and thickness $h/\sqrt{E}$ (solid line) of the steady supercritical axisymmetric intrusion as a function of the scaled radial distance $rE$ , given by (3.8). This solution illustrates both the radial thinning and the constant velocity obtained for $rE\gg 1$ . The constant far-field velocity is approached rapidly: $u/\sqrt{E}$ is within 4 % of its far-field value by $rE=2$ .

The time-dependent evolution of a radial intrusion, specified by the hyperbolic equations (3.1) and (3.2), requires two boundary conditions to be imposed at the source $r=r_{0}$ , because it is supercritical, since two characteristics leave this boundary (Whitham Reference Whitham1974). Both the flux condition (3.4) and the Bernoulli constant (3.7) used to determine the steady supercritical solution (3.8) are therefore specified at the source. We may identify the parameter $E$ in terms of a source Froude number, $\mathit{Fr}_{0}=2|u_{r}|/h$ at $r=r_{0}$ (where for this supercritical flow $\mathit{Fr}_{0}\geqslant 1$ ), from which we obtain

In the solutions that follow, the quantities $\mathit{Fr}_{0}$ and $r_{0}$ occur only through $E$ , and we may therefore specify without loss of generality either $\mathit{Fr}_{0}$ or $r_{0}$ . A natural choice for axisymmetric intrusions is to fix $\mathit{Fr}_{0}=1$ ; that is, to identify $r_{0}$ as the radius at which the flow is critical. In this case, from (3.11), we obtain $r_{0}=1/E$ .

In volcanic eruptions, the axisymmetric shallow intrusions that form are fed by a rising plume. Near the origin, the flow transitions from a near-vertical plume to a near-horizontal intrusion, and flow in this region is fully three-dimensional. In this section we have demonstrated that, once the intrusion has become shallow, the parameters at source that determine its behaviour are the volume flux and the parameter $E$ . While the volume flux is frequently calculated from observations of volcanic plumes, we note that $E$ may be more difficult to infer. In particular, our identification of $r_{0}=1/E$ as the radius at which the flow becomes critical is not the same as the commonly used ‘corner radius’ of a volcanic umbrella cloud (Sparks et al. Reference Sparks, Bursik, Carey, Gilbert, Glaze, Sigurdsson and Woods1997), a parameter characterising the radius at which the vertical flow in the rising plume transitions to a horizontal flow in the intrusion. The fully three-dimensional flow in the region at which the corner radius is defined means that there may be no simple relationship between the corner radius and the parameter $E=1/r_{0}$ which determines the behaviour of the intrusion. However, we will show that the growth rate of the intrusion at late times depends only very weakly on $E$ (the growth rate is proportional to $E^{1/8}$ ), indicating that accurate measurement of $E$ may not be necessary to obtain good predictions of the growth of volcanic intrusions.

3.2. Unsteady drag-free flow

We now examine the unsteady radial motion. We numerically integrate the unsteady radially symmetric equations (3.1), (3.2), with boundary conditions (3.4), (3.7) and (3.5). Here and throughout this paper we use the shock-capturing method of Kurganov & Tadmor (Reference Kurganov and Tadmor2000), which ensures that (2.4) is enforced at any discontinuities in the solution. We obtain the thickness field $h$ as a function of time (figure 2 a), along with the temporal evolution of the radial front $r_{f}(t)$ (figure 2 b). This evolution and the approach to steady state take a rather unusual form: rather than progressing as a straightforward similarity solution, as is common for many models of gravity current motion (see, for example, Ungarish Reference Ungarish2009), the time-dependent solution is steady everywhere except in a time-dependent region close to the front. The steady inner tail region is described exactly by the steady supercritical solution (3.8a,b ). The convergence of hyperbolic characteristics at the interface between this tail and the time-dependent frontal region means that the two regions are connected by an outward-propagating jump or shock (see, for example, Whitham Reference Whitham1974). We denote the radial position of this shock as $r_{s}(t)$ , which is determined as part of the solution (figure 2 a).

Figure 2.

Numerical time-dependent solutions of an axisymmetric intrusion in the absence of drag ( $\mathit{Fr}_{f}=1,E=1$ ). (a) The flow thickness $h$ , centred on the neutral buoyancy height $z=0$ , as a function of radial position at times $t=2,10,25,50,100$ . The radial positions of the front of the intrusion $r_{f}(t)$ and of the internal jump $r_{s}(t)$ are indicated for $t=100$ . (b) The radial position of the front of the intrusion $r_{f}$ and width of the time-dependent frontal region $r_{f}-r_{s}$ as a function of time. Dashed lines indicate the theoretical late-time predictions (3.24) and (3.25) for the current radius $r_{f}$ and the width of the time-dependent frontal region $r_{f}-r_{s}$ respectively.

Importantly, this two-region structure implies that the rate of radial spreading is rather different from what might have been anticipated on the basis of simple scaling arguments (see, for example, Lemckert & Imberger Reference Lemckert and Imberger1993; Woods & Kienle Reference Woods and Kienle1994). Such scaling arguments advance by requiring that the time dependence of the current length and flow variables are governed by similarity scalings. In a scaling argument, the linear increase in current volume due to the continuous source would imply that $r_{f}^{2}h\sim t$ , while at the intrusion front, $u_{r}\sim r_{f}/t$ is required for kinematic consistency and $u_{r}\sim h$ to match the imposed frontal Froude number. Together these would imply that $r_{f}\sim t^{2/3}$ . Our numerical solutions of the governing equations indicate that instead $r_{f}$ varies as a power of time that approaches $t^{0.75}$ at late times (figure 2 b), close to the rate $r_{f}\sim t^{0.72}$ found experimentally for large-scale continuously supplied intrusions by Lemckert & Imberger (Reference Lemckert and Imberger1993).

The form of the numerical solutions indicates that this discrepancy in growth rate between numerical solutions and scaling arguments is because the bulk of the motion away from the front does not scale in a manner determined by the conditions at the intrusion front, but is instead independent of time. The existence of time dependence only near the current front is entirely analogous to that found by Garvine (Reference Garvine1984) and Slim & Huppert (Reference Slim and Huppert2011) for axisymmetric gravity currents propagating over a horizontal boundary through a uniform environment due to a sustained source. In appendix A we show that the asymptotic arguments developed below may be applied also to such gravity currents. In fact the non-existence of the similarity solution for a sustained axisymmetric gravity current was established by Grundy & Rottman (Reference Grundy and Rottman1986), although they did not demonstrate what might occur instead. For both intrusions and gravity currents, we find numerically that the ordinary differential equations (ODEs) that result from assuming that the current is self-similar become singular in the interior of the flow, and cannot be integrated to satisfy the boundary conditions at both the source and the front of the intrusion.

We now explain why the similarity solution is not achieved and how this behaviour can be identified at the outset, without numerical solution of the governing equations. In contrast to the axisymmetric continuously supplied flows described above, gravity currents generated by a line source (e.g. Gratton & Vigo Reference Gratton and Vigo1994) do attain a similarity form, and the scalings of these can be obtained by simple scaling arguments. We therefore identify a fundamental difference between radial flows supplied by a point source (either gravity currents or intrusions) and ‘two-dimensional’ flows generated by a line source. This difference can be explained by dimensional analysis. If a unique time-dependent length scale can be constructed from time and the parameters that feature in the problem, a similarity solution for the current is expected (Barenblatt Reference Barenblatt1996). In this case there is exactly one dimensional group that can be constructed from the parameters of the problem. When more than one dimensional group can be constructed, a similarity solution may be obtained asymptotically at late times, but a non-similarity form is also possible.

The governing equations for an axisymmetric intrusion, (3.1) and (3.2), contain no parameters when redimensionalised and written in terms of $u_{r}^{\star }$ and $b^{\star }=N^{\star }h^{\star }$ : the parameters occur only in the boundary conditions. One dimensionless parameter $\mathit{Fr}_{f}$ enters into the problem in the boundary condition at the flow front (3.5), and at the source of an axisymmetric flow we write the two boundary conditions required for a supercritical flow as

These two parameters $[Q^{\star }N^{\star }]\equiv L^{3}T^{-2}$ and $[E^{\star }]\equiv L^{2}T^{-2}$ have differing dimensions, and consequently the only dimensionless parameter that can be constructed from them, $(Q^{\star }N^{\star })^{2}/(E^{\star 3}t^{\star 2})$ , is time-dependent. Dimensional analysis then indicates that the current radius at late times is given by $r_{f}=(Q^{\star }N^{\star }/E^{\star })f[\mathit{Fr}_{f},(Q^{\star }N^{\star })^{2}/(E^{\star 3}t^{\star 2})]$ , where $f$ is an arbitrary function of all the non-dimensional parameters in the problem. The occurrence of a time-dependent argument of the undetermined function $f$ means that the time dependence of the system is also undetermined, and a simple similarity solution is not necessarily expected.

In contrast, for an intrusion supplied by a line source of volume flux $Q_{l}^{\star }$ per unit width, the boundary conditions to be applied at source are

In this two-dimensional flow there is a single dimensional group specified by the parameters, which have dimensions $[Q_{l}^{\star }N^{\star }]\equiv [E^{\star }]\equiv L^{2}T^{-2}$ . We can then determine that the current length $x_{f}$ attains a similarity form $x_{f}=(Q_{l}^{\star }N^{\star })^{1/2}t^{\star }f(\mathit{Fr}_{f},Q_{l}^{\star }N^{\star }/E^{\star })$ , where $f$ is again a function of all the non-dimensional parameters in the problem, but in this case with no time-dependent parameters.

In contrast to the intrusions supplied by a source of constant flux, intrusions that result from an instantaneous release of a fixed volume of mixed fluid are governed by similarity solutions both in ‘two-dimensional’ and axisymmetric geometries (Ungarish Reference Ungarish2005; Ungarish & Zemach Reference Ungarish and Zemach2007). In such constant-volume releases, neither the dimensional governing equations (when written in terms of $u_{r}^{\star }$ and $b^{\star }=N^{\star }h^{\star }$ ) nor the boundary conditions involve a dimensional parameter: the only dimensional parameters occur in the specification of the initial conditions. In general, the parameters describing the initial conditions would not be expected to affect the current at late times. However, integrating the equation of conservation of mass across the current and applying the kinematic condition at the flow front implies that one of these parameters, the total current volume multiplied by $N^{\star }$ , is conserved throughout the time evolution of the current, and therefore plays a role in the buoyancy-inertial similarity solutions attained at late time. Importantly, no other parameter resulting from the initial conditions is conserved, so at late times the problem is described by only one dimensional parameter. In two-dimensional intrusions, the conserved quantity is $\int _{0}^{x_{f}^{\star }}b^{\star }\text{d}x^{\star }$ (where $x_{f}^{\star }$ is the location of the front of the current), which has dimensions $[L^{2}T^{-1}]$ , leading to a current that grows as $t^{1/2}$ (Ungarish Reference Ungarish2005). In axisymmetric intrusions, this conserved quantity is $\int _{0}^{r_{f}^{\star }}rb^{\star }\text{d}r^{\star }$ , with dimensions $[L^{3}T^{-1}]$ , leading to a current that grows as $t^{1/3}$ (Ungarish & Zemach Reference Ungarish and Zemach2007).

The results of a dimensional analysis for gravity currents flowing over a horizontal surface through an unstratified ambient are entirely analogous those obtained above for intrusions. As with intrusions, only a single dimensional parameter can be constructed for instantaneous releases of dense fluid in either two-dimensional or axisymmetric geometries, and for two-dimensional gravity currents supplied by a constant source of buoyancy, leading to similarity solutions in these three cases (Hoult Reference Hoult1972; Gratton & Vigo Reference Gratton and Vigo1994). For axisymmetric gravity currents resulting from a constant source of buoyancy, two dimensional parameters can be constructed and a non-similarity solution results (appendix A).

3.3. Details of the unsteady drag-free motion

We seek an analytical form for the non-similarity solution describing axisymmetric continuously supplied intrusions at late times that is suggested by dimensional analysis and by our numerical computations (figure 2). In constructing such an analytical form we obtain algebraically the observed $t^{3/4}$ growth rate of the current, and how this growth rate arises from the interaction of the two regions in the current, the time-dependent front and the steady inner tail. The jump at at $r=r_{s}(t)$ separating these two regions plays an important role in the dynamics of the motion. For this radial flow the jump conditions for mass and radial momentum (2.4) that act across this discontinuity become

where the shock speed $c=\text{d}r_{s}/\text{d}t$ and where $r_{s}^{+}$ and $r_{s}^{-}$ denote values of the radius $r$ immediately downstream and upstream of the shock respectively. On the source side of the discontinuity $(r<r_{s})$ the flow is described by the steady supercritical solution (3.8) and in particular when $r\gg r_{0}$ we find that

Since $h(r_{s}^{+})\gg h(r_{s}^{-})$ , we may deduce the leading-order contributions from the jump conditions (3.16), namely that $c=u_{r}(r_{s}^{+})$ and that $(1/12)h(r_{s}^{+})^{3}=u_{r}(r_{s}^{-})^{2}h(r_{s}^{-})$ .

The decrease in Froude number that occurs as fluid passes through the shock allows the steady supercritical flow, which has a Froude number that increases with increasing radial distance from source ( $\mathit{Fr}=2u_{r}/h\sim r$ at for $r\gg r_{0}$ ), to match onto the boundary condition of constant Froude number at the flow front (3.5). It is for this reason that only the supercritical steady solution (3.8) is relevant to the time-dependent flow. The Froude number in the subcritical steady solution (3.9) decreases with increasing radius ( $\mathit{Fr}\sim 1/r$ for $r\gg r_{0}$ ). If this steady subcritical solution were connected with a shock to the constant Froude number front condition (3.5), at large radial distances, fluid passing through this shock would move from a region in which $\mathit{Fr}\ll 1$ to a region in which $\mathit{Fr}$ is of order unity. Such an increase in Froude number is associated with production of energy at the shock (Whitham Reference Whitham1974), which is physically inadmissible. Thus, the supercritical steady solution is selected in time-dependent flows by the constant Froude number at the flow front (3.5).

We now seek the powers of time that govern the intrusion radius, shock position and flow variables in the front region. The momentum jump condition and asymptotic form of the steady inner tail solution (3.17) for $r\gg r_{0}$ imply a scaling on the flow thickness in the outer frontal region just downstream of the shock,

The dynamic and kinematic boundary conditions at the current front imply a scaling on the flow thickness at the front,

Assuming that the positions of the shock $r_{s}$ and intrusion front $r_{f}$ scale as the same power of time at late times ( $r_{f}\sim r_{s}$ ) and that the flow thickness likewise scales uniformly with time throughout the frontal region [ $h(r_{f})\sim h(r_{s})$ ] allows us to combine the scalings (3.18) and (3.19) to obtain that $r_{f}\sim r_{s}\sim t^{3/4}$ and that, in the frontal region, $u\sim h\sim t^{-1/4}$ . This scaling of $r_{f}(t)$ is consistent with the exponent determined by our numerical solutions (figure 2 b).

We therefore write $r_{f}=Kt^{3/4}+\cdots \,$ and $r_{s}=K_{s}t^{3/4}+\cdots \,$ , where $K$ and $K_{s}$ are constants to be determined, and pose similarity solutions for the thickness and velocity field within the frontal region of the form $h=t^{-1/4}H_{0}(r/r_{f})$ and $u_{r}=t^{-1/4}U_{0}(r/r_{f})$ . Two distinct current structures are possible with this ansatz. If $K_{s}<K$ , the width of the frontal region $r_{f}-r_{s}\sim t^{3/4}$ and the spatial variation of $h$ and $u$ within the frontal region is determined by the functions $H_{0}$ and $U_{0}$ . Alternatively, if $K_{s}=K$ , the width of the frontal region $r_{f}-r_{s}$ must grow more slowly than $t^{3/4}$ , and the thickness and velocity within the frontal region are spatially constant to leading order, taking the form $h=t^{-1/4}H_{0}(1)$ and $u=t^{-1/4}U_{0}(1)$ . We show now that this latter situation is the case.

The mass conservation equation (3.1), when integrated across the frontal region, yields

At leading order we have found that $c=u(r_{s}^{+})$ , but by substituting the leading-order expressions and asymptotic form of the steady solution (3.17) into the jump condition for mass (3.16a ), we obtain that at the next order,

Thus, on integrating (3.20), we obtain

Using the leading-order expressions for $r_{s}$ , $r_{f}$ and $h$ in (3.22), we find that the powers of time on either side of the equation can balance only if $K_{s}=K$ .

Thus, to leading order the shock and the front move at the same speed. Substituting the leading-order expressions $h=t^{-1/4}H_{0}(1)$ and $u=t^{-1/4}U_{0}(1)$ into the kinematic and dynamic boundary conditions at $r=r_{f}$ , we obtain $U_{0}(1)=3K/4$ and $H_{0}(1)=3K/(2\mathit{Fr}_{f})$ respectively. From the momentum jump condition (3.16b ),

which determines that the coefficient multiplying the intrusion radius $r_{f}=Kt^{3/4}$ is

The width of the frontal region at late times can be determined by substituting the similarity forms into (3.22), from which we obtain

This predicted width, and the current radius (3.24), compare well with those calculated from our time-dependent numerical computation (figure 2 b). We note that at late times, (3.22) implies that the volume of the fluid within the outer frontal region of the intrusion scales as $r_{f}(r_{f}-r_{s})h\sim t$ , whereas the volume of fluid within the tail region $\int _{0}^{r_{s}}rh\,\text{d}r\sim t^{3/4}$ . Thus, in the asymptotic regime of late time, the bulk of the intrusion volume is contained within the time-dependent frontal region (cf. figure 2).

We demonstrate in appendix B that the long-time behaviour of the current can also be derived by posing an asymptotic expansion for the flow variables. In doing so we obtain a more accurate description of the flow at long times, by calculating next-order corrections to the leading-order solution above.

3.4. Unsteady axisymmetric flow with drag

Our solutions so far have neglected the effect of drag on the radial intrusions. The neglect of drag is an appropriate dynamical regime provided the drag force $(C_{D}u_{r}|u_{r}|)$ is much smaller than the inertia $\partial (u_{r}h)/\partial t$ . By comparing the magnitudes of these terms in the frontal region at late times, we deduce that drag may be neglected here if $C_{D}\mathit{Fr}_{f}t\ll 1$ . Numerical time-dependent solutions of the axisymmetric equations with drag (figure 3 a) are consistent with this view, with solutions with $C_{D}\mathit{Fr}_{f}t\ll 1$ closely resembling those in the drag-free case (figure 2 a).

Figure 3.

Time-dependent solutions of an axisymmetric intrusion including effects of drag. (a) Illustration of the flow thickness $h$ centred on the neutral buoyancy height $z=0$ , with $C_{D}=0.01$ and parameters otherwise as in figure 2 at times $t=2,10,25,50,100,200$ . (b) Current front position $r_{f}$ as a function of time, for $C_{D}=0.1$ , 0.01 and 0.001. The dashed line indicates the theoretical growth rate at late times of the current with $C_{D}=0.1$ . (c) Numerical solutions of the time-dependent equations with $C_{D}=0.1$ at time $t=2\times 10^{3}$ , rescaled to similarity variables (solid curves), compared with solutions to the similarity ODEs (3.28) for the drag-dominated regime (dashed curves).

In dimensional form, we therefore have that drag becomes significant when $t^{\star }\gg 1/(N^{\star }C_{D}\mathit{Fr}_{f})$ . This time is independent of the volume flux of the eruption. The value of $C_{D}$ for an intrusion at high Reynolds number is poorly constrained, but by analogy with the turbulent drag that results from entrainment (e.g. Fernando Reference Fernando1991) we anticipate values between 0.001 and 0.1. With $C_{D}=0.01$ the estimate above, with typical atmospheric stratification of $N^{\star }=0.01~\text{s}$ , indicates that drag may become important after approximately 3 h.

When $C_{D}\mathit{Fr}_{f}t\approx 1$ (for example, the solutions at $t=50$ , 100 and 200 in figure 3 a) the effects of drag become significant. In figure 3(a) the reduction in velocity in the steady supercritical tail due to drag results in a thicker intrusion in this region, compared with the drag-free case illustrated in figure 2(a). The time-dependent frontal region is both wider and thinner than in the drag-free case, resulting in an intrusion with a much smaller change in thickness across the jump linking the two regions. The width of the frontal region grows rapidly with time when $C_{D}\mathit{Fr}_{f}t>1$ , as drag becomes increasingly significant, with this frontal region at late times forming the majority of the intrusion.

When $C_{D}\mathit{Fr}_{f}t\gg 1$ a new dominant momentum balance is established between the radial pressure gradient and the drag. Then, similarly to gravity current motion in a uniform environment (Hatcher, Hogg & Woods Reference Hatcher, Hogg and Woods2000; Hogg & Woods Reference Hogg and Woods2001), the flow undergoes a transition so that the dimensionless governing equations for the radial motion are

These parabolic equations are subject to the source condition (3.4). In this drag-dominated regime, there is a new class of similarity solutions for the ensuing motion, and scaling can be usefully employed to determine the similarity exponents. To this end, we note from integrating (3.26a ) across the intrusion that mass conservation demands $r_{f}^{2}h\sim t$ , while from (3.26b ) the dynamical balance requires $h^{3}/r_{f}\sim C_{D}r_{f}^{2}/t^{2}$ . Thus, we deduce that the similarity variable $r_{f}\sim C_{D}^{-1/9}t^{5/9}$ and seek solutions of the form

where ${\it\eta}=r/r_{f}(t)$ and ${\it\kappa}$ is a dimensionless constant to be determined. On substitution of (3.27) into the governing equations (3.26), we obtain

where the prime denotes differentiation with respect to ${\it\eta}$ . These are subject to boundary conditions $\mathscr{U}(1)=5/9$ and $\mathscr{H}(1)=0$ , where the former imposes kinematic consistency at the front and the latter is the Froude number condition in the regime $C_{D}t\gg 1$ . Finally, the specification of the source flux is most easily enforced by integrating the mass conservation equation (3.26a ) from $r=r_{0}$ to $r=r_{f}$ , which yields

because $r_{0}/r_{f}(t)\ll 1$ when $C_{D}t\gg 1$ . Only one boundary condition at source is required for these subcritical drag-dominated flows. The ordinary differential equations (3.28) are readily integrated numerically (for example, using a Runge–Kutta method), and applying (3.29) we find that ${\it\kappa}=1.222$ . The profiles of $\mathscr{H}$ and $\mathscr{U}$ , illustrated in figure 3(c), closely match the numerical solution to the time-dependent equations (3.1) and (3.2) at late times. In the similarity regime both the flow thickness and the velocity decay monotonically with increasing radius, with the intrusion exhibiting a rounded flow front.

The existence of a similarity solution in this drag-dominated regime (in contrast to the non-similarity solution obtained in the drag-free regime) may be anticipated from dimensional analysis. By writing (3.26a,b ) in terms of a new variable $hN^{2/3}$ , we find that the problem can be written in terms of only a single dimensional parameter, $N^{\star 2}Q^{\star 3}$ , with dimensions $[N^{\star 2}Q^{\star 3}]\equiv L^{9}T^{-5}$ . Thus, as in the analysis for an intrusion generated by a line source in § 3.2, the only length scale is a similarity length scale, which in this case grows as $t^{5/9}$ , and a similarity solution to the equations must result.

We note that for ${\it\eta}\ll 1$ the solution has the form $\mathscr{H}({\it\eta})\sim {\it\eta}^{-1/5}+\cdots \,$ and $\mathscr{U}({\it\eta})={\it\eta}^{-4/5}+\cdots \,$ , and thus the original dimensionless variables $h$ and $u_{r}$ are independent of time for $r\ll r_{f}$ . To first order in $r/r_{f}$ , mass conservation (3.26a ) then implies that $ru_{r}h$ is constant, from which we find

These leading-order terms were obtained by Baines (Reference Baines2013) under the assumption of a steady flow. We have shown that that once drag has become dominant, this steady solution is in fact attained by the attracting time-dependent similarity solution, but only in the regime $r\ll r_{f}$ . Thus, the flow is time-dependent close to the front of the intrusion, but a steady intrusion is ‘laid down’ asymptotically, closer to the source.

4.1. Downwind spreading regimes

Our calculation of current width spreading as $w\sim x^{1/3}$ far downstream differs from the prediction of $w\sim x^{1/2}$ of Bursik et al. (Reference Bursik, Carey and Sparks1992a ). Their prediction is made under the assumption that the downstream speed of the intrusion is equal to the ambient wind speed $U$ due to drag, but that the lateral spreading far downstream is in an inertia-buoyancy regime, governed by the constant Froude number condition at the flow extremities, rather than being governed by the balance between lateral pressure gradients (4.3) and drag. Drag becomes increasingly important in this lateral spreading at late times (or, equivalently, far downwind of the source), and thus we expect a transition some distance downwind from an intrusion spreading in the inertia-buoyancy regime of Bursik et al. (Reference Bursik, Carey and Sparks1992a ) to one spreading in the drag-dominated regime, as described in § 4. We seek the location of this transition.

The width of a wind-dominated intrusion in the inertia-buoyancy spreading regime is given by Bursik et al. (Reference Bursik, Carey and Sparks1992a ), up to an undetermined parameter, postulated to be of order unity, depending on the flow geometry. We now calculate this parameter within the context of the model proposed by Bursik et al. (Reference Bursik, Carey and Sparks1992a ) in which the downstream intrusion velocity is exactly equal to the wind speed, $u=U$ . Under this assumption the steady governing equations in the absence of drag become

and

Identifying the quantity $U\partial /\partial x$ with the temporal derivative in a frame moving downstream at velocity $U$ , as Bursik et al. (Reference Bursik, Carey and Sparks1992a ) do, we see that (4.8) and (4.10) are simply the equations of mass and lateral momentum conservation for a one-dimensional intrusion, with the dynamic boundary condition $v=\mathit{Fr}_{f}h$ at the flow front. The width of the steady laterally spreading intrusion can therefore be determined by solutions to the unsteady one-dimensional system (4.8) and (4.10), which has been studied previously in the context of intrusions resulting from the instantaneous release of a constant volume of mixed fluid. At late times a similarity solution describes the spreading (Ungarish Reference Ungarish2009), and predicts a thickness profile for the intrusion

where the intrusion half-width $w_{i}(x)$ is given by

and $\tilde{C}$ is a function of the frontal Froude number. We note that under the scalings implied far downstream by (4.11) and (4.12), the term on the right-hand side of (4.9) is asymptotically small, and thus this equation is also satisfied by any solution of (4.8).

Integrating the thickness profile (4.11) across the width of the flow, and requiring, as in (4.4), that the flux through this profile is equal to the source flux gives

Evaluating this integral for the value of the Froude number suggested by Ungarish (Reference Ungarish2009), $\mathit{Fr}_{f}=1.19/\sqrt{2}$ , we then obtain $\tilde{C}=1.926$ . It is notable that the thickness profile in the buoyancy-inertial spreading regime (4.11) is qualitatively very different from that found in the drag-dominated regime (4.6), with the minimum intrusion thickness occurring along the centreline, $y=0$ .

Equating the half-width $w_{i}$ given by (4.12) with the half-width predicted for the drag-dominated current (4.6), we find that the spreading becomes drag dominated at an approximate distance

downwind of the intrusion source. As we have noted, the value of $C_{D}$ for an intrusion at high Reynolds number is very poorly constrained, but using $C_{D}=0.01$ with realistic atmospheric values of $N^{\star }=0.01~\text{s}^{-1}$ and $U^{\star }=10~\text{m s}^{-1}$ , (4.14) indicates that the transition from inertia-limited to drag-limited spreading occurs in the region of 200 km downstream of the source. While the uncertainty in $C_{D}$ limits the accuracy of this estimate, the order of magnitude suggests that both regimes are likely to be observable in volcanic plumes, which can reach over 1000 km downwind (Pouget et al. Reference Pouget, Bursik, Webley, Dehn and Pavolonis2013). It is notable that the flux $Q^{\star }$ does not enter into the redimensionalised form of (4.14), $x_{t}^{\star }\approx 2.0U^{\star }/(C_{D}N^{\star })$ ; this estimate is therefore independent of the plume flux.