A dam break driven by a moving source: a simple model for a powder snow avalanche

We study the two-dimensional, irrotational flow of an inviscid, incompressible fluid injected from a line source moving at constant speed along a horizontal boundary, into a second, immiscible, inviscid fluid of lower density. A semi-infinite, horizontal layer sustained by the moving source has previously been studied as a simple model for a powder snow avalanche, an example of an eruption current, Carroll et al. (Phys. Fluids, vol. 24, 2012, 066603). We show that with fluids of unequal densities, in a frame of reference moving with the source, no steady solution exists, and formulate an initial/boundary value problem that allows us to study the evolution of the flow. After considering the limit of small density difference, we study the fully nonlinear initial/boundary value problem and find that the flow at the head of the layer is effectively a dam break for the initial conditions that we have used. We study the dynamics of this in detail for small times using the method of matched asymptotic expansions. Finally, we solve the fully nonlinear free boundary problem numerically using an adaptive vortex blob method, after regularising the flow by modifying the initial interface to include a thin layer of the denser fluid that extends to infinity ahead of the source. We find that the disturbance of the interface in the linear theory develops into a dispersive shock in the fully nonlinear initial/boundary value problem, which then overturns. For sufficiently large Richardson number, overturning can also occur at the head of the layer.

We study the two-dimensional, irrotational flow of an inviscid, incompressible fluid injected from a line source moving at constant speed along a horizontal boundary, into a second, immiscible, inviscid fluid of lower density. A semi-infinite, horizontal layer sustained by the moving source has previously been studied as a simple model for a powder snow avalanche, an example of an eruption current, Carroll et al. (Phys. Fluids, vol. 24, 2012, 066603). We show that with fluids of unequal densities, in a frame of reference moving with the source, no steady solution exists, and formulate an initial/boundary value problem that allows us to study the evolution of the flow. After considering the limit of small density difference, we study the fully nonlinear initial/boundary value problem and find that the flow at the head of the layer is effectively a dam break for the initial conditions that we have used. We study the dynamics of this in detail for small times using the method of matched asymptotic expansions. Finally, we solve the fully nonlinear free boundary problem numerically using an adaptive vortex blob method, after regularising the flow by modifying the initial interface to include a thin layer of the denser fluid that extends to infinity ahead of the source. We find that the disturbance of the interface in the linear theory develops into a dispersive shock in the fully nonlinear initial/boundary value problem, which then overturns. For sufficiently large Richardson number, overturning can also occur at the head of the layer.

Introduction
A powder snow avalanche is a flowing, fluidised suspension of fine snow, fed by new snow that is ejected from the underlying surface due to large pressure gradients near the head of the avalanche, sometimes referred to as an eruption current (Sovilla, Burlando & Bartelt 2006;Louge, Carroll & Turnbull 2011). This is in contrast to more commonly studied examples of gravity currents, which are formed by the release of a fixed mass of the particulate substance that then moves under the action of gravity without further replenishment (Simpson 1997). The presence of a source of new particulate material close to the head of the flow fundamentally changes its dynamics, and may allow an eruption current to propagate further than a simple A unit source at x = y = 0 in a unit flow from right to left when the two fluids have equal density, given by the velocity potential φ = −x + (1/π) log r and corresponding streamfunction ψ = −y + θ/π. The streamline that meets the stagnation point on the boundary at x = 1/π, y = 0 and which represents the position of the free surface, x = y cot πy (equivalently y = θ/π), is shown as a broken line.
gravity current. In Carroll, Turnbull & Louge (2012), the authors consider a sequence of two-dimensional models of powder snow avalanches, in each of which the snow that erupts into the current is modelled as either a point source in an inviscid, irrotational flow or as a numerical approximation of a point source within a large eddy simulation of the variable density Navier-Stokes equations. In this paper we revisit the first of these models, in which the fluidised snow is modelled as inviscid and incompressible, and is assumed to flow irrotationally from a moving point source surrounded by stationary, ambient, immiscible, inviscid fluid of lower density. When the fluid densities are equal, in a frame of reference moving at velocity U > 0 with the point source of strength q at its origin, the flow is given in polar coordinates (r, θ ) by the velocity potential φ = −Ur cos θ + q π log r, (1.1) in 0 θ π. Some typical streamlines are shown in figure 1, and are equivalent to those around a Rankine body given by in Cartesian coordinates (x, y) (see, for example, Batchelor 1967). In terms of the problem that we are interested in, equation (1.2) gives the dividing streamline between fluid 2 and fluid 1, which originates at the source, when the fluids have equal density. We find that the flow in this nonlinear, two fluid, free boundary problem is more complex than is suggested by the results presented in Carroll et al. (2012). Unless the fluids have equal density, no steady solution exists because of the likely nonexistence of a steady local solution in the neighbourhood of the contact point (at y = 0). We therefore analyse the dynamics of the flow by studying an initial value problem in which the initial velocity field is given by (1.1) and the fluids are initially separated by the free surface given by (1.2). Later in the paper, we also consider the flow when the initial position of the free surface is a nearby streamline of this flow, for reasons that will be become clear as we progress. The difference in density between the fluids drives the subsequent development of the flow and motion of the free surface. Although it may seem more natural to study what happens when a source of slightly denser fluid is injected ab initio, we are interested in how the dynamics of this flow is related to fully developed powder snow avalanches, for which the initial conditions that we have chosen are more relevant.
In § 2, we set up the governing equations for this initial/boundary value problem. We study the problem in the limit of small density difference in § 3, and show that the head of the layer of fluid 1 (to the right of the source) swells and the contact point moves rapidly forwards close to the initial stagnation point. In addition, a disturbance propagates downstream at the speed of the source, which is a stationary initial transient in a frame of reference fixed in the boundary. In § 4 we examine the initial behaviour of the flow, in particular near the stagnation point where the interface meets the solid surface. We find that this is locally a dam break, which means that we need to use the method of matched asymptotic expansions, in the same manner as King & Needham (1994) and Uddin & Needham (2015), to determine the initial dynamics of the flow. Close to the contact point, we find that the free boundary takes a spatially oscillatory form, with a wavelength that tends to zero as y → 0. We also find that in this region the Boussinesq approximation,ρ 1 andRi = O(1), which it is natural to consider in the outer flow (although it is not relevant to avalanches), is a singular limit. In § 5, we discuss numerical solutions of the fully nonlinear free boundary problem, computed using an adaptive vortex blob method, which allows us to study the nonlinear development of the initial disturbance. We find that nonlinear effects (Kelvin-Helmholtz instability) can cause the disturbance that propagates along the layer to overturn. For non-zero Richardson number, we regularise the flow close to the contact point by using an initial free surface with a thin layer of fluid 1 that extends to infinity ahead of the source. In this case, overturning can also occur at the head of the layer for sufficiently large Richardson number.

Governing equations
Using the source velocity, U, and the flux from the source, q, we can form length and time scales, L ≡ q/U and T ≡ q/U 2 . We use these to non-dimensionalise the governing equations in a frame of reference moving with the source. The initial/boundary value problem that we will study, defined by (2.2)-(2.8) below, contains just two dimensionless parameters constructed from these quantities: the dimensionless density difference and the Richardson number, In a powder snow avalanche, the density of the fluidised snow is typically between two and ten times greater than that of the surrounding air, the height of the avalanche is between 20 and 50 m, and typical propagation speeds lie between 15 and 50 m s −1 . This leads to the approximate parameter ranges 1 ρ 10, 0.1 Ri 20 and 0.1 R i 2. The Richardson number quantifies the importance of the effective gravitational force on the flow compared to that of the relative acceleration of the two fluids (the Bernoulli effect).
The two-dimensional domain of solution is y 0, −∞ < x < ∞ in Cartesian coordinates. We denote by D i (t) the region occupied by the fluid of density ρ i with velocity potential φ i (x, t), where t 0 is time, x = (x, y) and i = 1 or 2. The boundary ∂D(t) between these two regions is given by x = X(s, t) ≡ (X(s, t), Y(s, t)), where s is the arclength along the free surface, with s = 0 at the point x = X c (t) where ∂D meets the boundary, y = 0. The irrotational flow of the two incompressible, inviscid fluids is governed by with the kinematic boundary condition n · ∂X ∂t = n · ∇φ 1 = n · ∇φ 2 at x = X(s, t), where n is the outward unit normal to D 1 , along with the dynamic boundary condition, Note that we have set the Bernoulli constant to zero, but that this is somewhat arbitrary, since the interface does not extend to infinity upstream. The no penetration condition at the solid boundary is Although (2.3) shows that the normal derivatives of the potentials, and hence normal velocities, are continuous at the boundary, the potentials and their tangential derivatives (the tangential velocities) need not be continuous in this inviscid model. We must also include the source of fluid 1 at the origin, φ 1 ∼ 1 π log r as r ≡ x 2 + y 2 → 0, (2.6) and the upstream flow of fluid 2 in this frame of reference, Finally, we take the initial position of the free surface to be the Rankine surface for the flow with equal densities, given by (1.2), so that X = Y cot πY when t = 0, (2.8) and hence when t = 0, D 1 is the region 0 y Y for x 1/π, with D 2 the remainder of the upper half-plane. Although this is not an entirely obvious choice of initial condition -the fluids do not have equal density, so this is not an equilibrium solution, and the larger the density difference, the more rapidly the interface between the fluids evolves away from the initial condition -it has several advantages over other, more arbitrary choices, including: (i) it consists of a semi-infinite layer of the denser fluid, a geometry that we must preserve in order to be relevant to powder snow avalanches, (ii) it is close to an equilibrium solution when the density difference is small, so we can use a weakly nonlinear analysis to obtain insight into the development of The position of the free surface forρ = 0.1 andRi = 0, 1, 5 and 10 when t = 1, calculated from the asymptotic solution forρ 1. The initial position of the free surface is shown as a broken line.
the flow, (iii) it is easily modified to include a thin layer of the denser fluid stretching to infinity ahead of the source (see § 5), (iv) it is well defined and easy to implement both analytically and numerically. Although we have not tried to use any other initial conditions, it seems reasonable that the results below capture, at least qualitatively, the nature of the flow that would develop from similar initial conditions. Note also that in § § 3 and 4, we will usually assume (correctly) that the free surface can be written as a single-valued function in polar coordinates, r = R(θ , t), with corresponding initial condition R(θ , 0) = θ/π sin θ , or as y = Y(x, t) with x = Y(x, 0) cot πY(x, 0) or as x = X(y, t), with X(y, 0) = y cot πy, in Cartesian coordinates. In § 5, we will use a Lagrangian representation of the free surface, which can accommodate overturning, as part of our numerical solution method for the full initial/boundary value problem. Although φ 1 = φ 2 when t = 0, as the inviscid, irrotational flow develops, a discontinuity in the tangential derivatives of these potentials, and hence a tangential velocity discontinuity, develops at the free surface. Locally, this will resemble planar, inviscid, two fluid flow with a velocity discontinuity. Flows like this are known to be subject to Kelvin-Helmholtz instability (for example, Batchelor 1967), which usually leads to the development of a curvature singularity in finite time (for example, Moore 1978;Cowley, Baker & Tanveer 1999), which can be regularised by adding some extra physics to the model, usually surface tension, although this is not the physics that we will add in this paper. We will address this issue in detail in § 5.
3. Asymptotic solution for flows with small density difference (ρ 1,Ri = O(1)the Boussinesq approximation) Although we have seen thatρ is not small for powder snow avalanches, it is natural to consider the Boussinesq limit,ρ 1 withRi = O(1), as a way to gain insight into the dynamics of the flow and provide a test for numerical solutions of the full problem. As we shall see, in this limit the head of the layer swells and a depression propagates back along the layer at unit speed. We will see in § 5 that this is also how the flow behaves whenρ is not small. We will also discover that there is a singularity in the shape of the free surface close to the contact point (a dam break), which persists when ρ is not small. The mathematical details of our analysis are given in appendix A. Figure 2 shows the position of the free boundary when t = 1 forρ = 0.1 and various values ofRi, calculated using the asymptotic solution forρ 1. This somewhat large value ofρ has been chosen to allow the position of the interface to be seen more clearly. In each case, the layer swells at the head, with the size of the swelling increasing withRi andρ, whilst a depression propagates back down the layer at unit speed. Note that whenRi = 0 (the zero gravity case), the free surface remains perpendicular to the boundary, y = 0, and the speed of the contact point asymptotes to zero as t → ∞. WhenRi > 0, a thin layer of fluid 1 protrudes into fluid 2 close to the boundary, with the shape of the free surface very different to that in the zero gravity case. This is because there is a logarithmic singularity in the position of the free surface as y → 0, which we have not dealt with explicitly in our evaluation of the integral in (A 18), and is therefore removed by rounding error when y = 0. The behaviour close to the boundary when Ri = O(1) and t 1 is the subject of § 4, and includes the limitρ 1, which we have studied in this section, as a subcase. In order to see why there is a singularity, consider the dynamic boundary condition, equation (2.4). Initially, the free surface is perpendicular to the boundary, y = 0, so we can differentiate (2.4) with respect to y, and use the no penetration condition, equation (2.5), to show that the left hand side is then zero. The right-hand side is, however, non-zero, and it is this incompatibility in the boundary conditions that leads to a singularity in the outer flow when t 1. This incompatibility is reconciled in a small inner region close to the point of contact of the free surface with the boundary which, crucially, moves with an algebraically small velocity and a logarithmically large acceleration as t → 0. In § 5 we determine an equation for the motion of the contact point, equation (5.7), that holds for all t > 0. Figure 3 shows the free surface at later times for various values ofRi andρ, again calculated using the asymptotic solution forρ 1. The depression of the free boundary can be seen to continue to propagate at unit speed whilst deepening, steepening and eventually overturning (although the position of the free boundary remains a single-valued function r = R 0 (θ) +ρR(θ , t)). The larger the value ofRi, the deeper the depression (note the different choices ofρ in each panel of figure 3). As the depression deepens and the slope of the boundary increases, this linear theory breaks down. As we shall see in § 5, nonlinearity leads to a dispersive shock at the rear of the depression (see figure 10), and the eventual formation of a singularity in the curvature. Ahead of the depression, the free boundary asymptotes to a locally steady solution with thickness larger than that of the original layer. Note however that, as discussed above, forRi > 0, the position of the free surface is logarithmically singular as y → 0. The bounded solution computed in this outer region is an artefact of the method that we have used to approximate the integral in (A 18), with rounding error cutting off the singularity. This indicates how easy it is to compute spurious numerical solutions for this problem, since the logarithmic singularity is so weak. The solution when t 1 withρ = O(1) and Ri = O(1) consists of two asymptotic regions: an outer region, where the solution is a small perturbation of the initial conditions, and an inner region, close to the point of contact of the free surface with the boundary, y = 0, where the logarithmic singularity in the outer flow is resolved. In each region, it is most convenient to represent the free surface as x = X(y, t). Note that for this problem, we find that the Boussinesq approximation (ρ = 0 with Ri > 0) is a singular perturbation. The acceleration of the two fluids in the inner region is crucial in resolving the outer singularity, and if the two fluids are treated as having equal density in the inertial terms in (2.3), this acceleration is not associated with a restoring force. As noted earlier, the flow in the inner region is effectively a two-fluid dam break, and we use the methods developed by King & Needham (1994) and Uddin & Needham (2015) to determine the asymptotic solution as t → 0. In addition to the solution in each of two asymptotic regions, constructed and described in appendix B, it is straightforward to construct a composite solution, and figure 4 shows a typical example. Note the huge difference in scales between the main and inset figures, even for the relatively large value, t = 0.1, which illustrates the complex structure of the asymptotic solution. The displacement of the contact point is of O(t 2 (−log t)) as t → 0, and this singular initial acceleration gives rise to oscillations in the free surface that are of decreasing amplitude and wavelength such that the slope is bounded but the curvature is singular as y → 0. See appendix B for a more detailed discussion.

Regularisation by a precursor layer
The natural next step in this investigation is to solve the full, nonlinear initial/boundary value problem numerically. We should, however, pause to take stock of what we now know about the behaviour of the solution. We have discovered that, when Ri > 0, the initial flow near the head of the fluid layer is highly spatially oscillatory and develops in a region with length scale of O(t 2 ) as t → 0, with a singular curvature at the contact point. The intricate, multiscale nature of the free surface close to the boundary when t 1 (see figure 4) means that it is unclear whether it is possible to start a simulation using the asymptotic solution as the initial condition. We did attempt to do this, but were unable to develop a numerically stable solution method. In addition, since the asymptotic solution discussed in § 4 is the leading-order term in an expansion in powers of (−log t) −1 , its accuracy at times small, and yet large enough to use to initialise a numerical simulation, is low.
In the similar, single fluid dam-break flow studied in Uddin & Needham (2015), the inclusion of weak surface tension regularised the solution at small times. This was also the case in the problem studied in Billingham et al. (2017), where a numerical solution showed that the oscillatory small time solution emerges after a short, surface-tension-dominated transient. Whilst we could do the same here, including surface tension is not in the spirit of a simplified problem that is meant to be related to powder snow avalanches, for which the small scale behaviour close to the free boundary is better characterised as particle-laden, multiphase turbulence. We therefore decided, in addition to the regularisation provided by the vortex blob method (see discussion below), to modify the initial condition so that a thin layer of fluid 1 (stationary as x → ∞ in a frame of reference fixed in the solid boundary) lies ahead of the main layer. This is easily achieved by using the streamline given by as the initial position of the free boundary. In this way, the single fluid solution remains correct whenρ = Ri = 0, and has Y → as X → ∞. Figure 5 shows some numerical solutions for the typical case Ri = 5,ρ = 2, calculated for a sequence of successively smaller values of , in a frame of reference fixed in the boundary, so that the source (shown as a circle) moves to the right at unit speed (see § 5.3 for details of the numerical method used). The features displayed in these simulations, namely a depression of the free surface near the original position of the source, overturning of the upper surface of the layer, and overturning at the head of the layer, are, as we shall see, typical. Note that the stopping criterion was that the number of marker points required for the resolution of the free surface that we used reached 1500, which is a good indicator that the free surface is about to pinch off. Figure 5 provides strong evidence that the value = 10 −3 , which we have used in the simulations shown below, is a reasonable choice, as the structure of the flow and the time at which the free surface pinches off appears no longer to be changing as decreases past this value. We have not attempted to construct the asymptotic solution in the limit → 0 in this paper. Since this regularisation removes the contact point, we can also consider looking for steady solutions. Now Y → and φ i ∼ −x as x → ∞, so we must modify the Bernoulli constant in (2.4) to give We tried to find steady solutions using a numerical method based on the steady version of that presented below in § 5.3 and continuation from the known solution whenρ = Ri = 0. We tried continuation inρ, and also in an artificial parameter multiplying the source term. We were, however, unable to obtain converged solutions. If we investigate this further by considering the asymptotic solution forρ 1, the equivalents of (A 6) and (A 19) show that now This leads to a non-integrable singularity in the vortex sheet representation, equation (A 18), which shows that a finite solution does not exist, consistent with our failed attempts to use numerical continuation starting fromρ = 0. All of this strongly suggests that no steady solution of the regularised problem with > 0 exists, and hence that the steady solution φ i = −x + (1/π) log r is structurally unstable to small perturbations in the density of the fluid below a streamline. This does not, however, mean that the unsteady solution that evolves from initial conditions based on this special solution are not interesting and of possible relevance to powder snow avalanches.

The adaptive vortex blob method
It is well known that periodic, two-fluid, inviscid, irrotational flow with a vortex sheet is ill posed (Ebin 1988) and that Moore singularities develop in finite time (Moore 1978). When surface tension is included, the singularity is avoided, and instead the free surface forms Kelvin-Helmholtz rolls. Spectrally accurate simulations of this roll-up have been performed on geometrically simple problems with periodic initial conditions, regularised either by surface tension (for example, Hou, Lowengrub & Shelley 1994) or by using a vortex blob method (for example, Baker & Pham 2006). We therefore choose to use the latter method, in which the vortex sheet is replaced with a distribution of vorticity.
There are some remarks that we should make that arise from the aperiodic nature of the flow that we are interested in. Firstly, we are not concerned here with the details of the Kelvin-Helmholtz roll-up process. We want to understand the bulk dynamics of the flow, and determine where roll-up starts to occur, but appreciate that the details of the roll-up will depend on the regularisation used. Secondly, in previous work using the vortex blob method, the grid spacing in the periodic domain of solution was constant. We will use an adaptive boundary integral method, with a time-dependent, non-uniform grid size that allows us to resolve regions of the boundary where small scale features develop. We must therefore decide how to choose L δ , the length scale over which the vorticity is distributed. We will consider two approaches. In Method 1, we take the length scale to be the grid size, L δ = s. This is equivalent to a choice taken by Baker & Beale (2004), which was shown to lead to accurate results for periodic flow without adaptive regridding. In Method 2, we take L δ = max( s, δ), with δ > 0 a constant. Here, δ is a regularisation parameter that gives the length scale over which vorticity is smeared, which we can think of as a physical quantity that characterises sub-element scale processes, such as turbulence and particle-particle interactions. Using Method 1 (which we note is equivalent to Method 2 with δ = 0), we are able to compute the development of a Moore singularity up to large curvatures, but, because we adaptively regrid the boundary to resolve the local curvature, the singularity is not removed, and we cannot compute beyond it. The solution is plotted in a frame of reference where the outer fluid is initially stationary at infinity and the source (indicated as a small circle) moves to the right at unit speed. Note the different scales for the axes in each panel. The simulation with δ = 0 was terminated when the maximum curvature was greater than 1000, and for δ > 0, just before pinch off of the interface.
The effect of varying δ in the typical example, Ri = 5,ρ = 2, is shown in figure 6 (note the different scales on the axes). When δ = 0, the solution terminates at a developing singularity in the curvature near the head of the layer, as shown in figure 7, where we compute until the curvature reaches 1000. In contrast, using Method 2, the Moore singularity is regularised, so we can integrate beyond its formation time and compute an interface that rolls up and overturns, plotted in figure 6 just before the free surface pinches off. These rolls start to form at the same time and place as in the solution for δ = 0. As δ increases, the onset of overturning is delayed and the length scale of the rolls increases. In the simulations shown below, we have used δ = 0.05.

Numerical method
In order to solve the initial/boundary value problem given by (2.2)-(2.8), it is natural to use a boundary integral method. Consider a set of marker points in the boundary, ∂D, parameterised by α, that move with velocity U ≡ Ut + (u · n)n, so that where d/dt is the Lagrangian time derivative. The normal component of this velocity matches that of the two fluids, but the tangential component is determined, as we shall see below, by specifying the spacing of the surface marker points. We define the unnormalised vortex sheet strength to be γ ≡ µ α (α, t). If γ is known, the velocity of the fluid is given, using derivatives of the Plemelj formula (B 13), and an image of the free surface in the solid boundary, by u, v). The first two terms on the right-hand side of (5.5) are the unit flow and unit source, which means that γ = 0 when t = 0. Following Baker, Meiron & Orszag (1982) and Hou et al. (1994), the dynamic boundary condition can be manipulated to find the evolution equation for γ , as where s(α, t) is arclength and W is the fluid velocity at the interface given by the principal value integral (5.5). Note that n · W is continuous at the boundary, but t · W is not, and jumps by an amount γ /s α across the boundary. Note also that, since the first step in the lengthy derivation of (5.6) is differentiation of (5.2) with respect to α, the Bernoulli constant does not appear. The Bernoulli constant in (5.2) simply reflects the upstream boundary conditions, which are satisfied by the initial conditions in the evolution problem that we solve here. In all the solutions shown below, γ → 0 as α → ±∞, as we would expect since the far field is in equilibrium both up-and downstream.
As an aside, this form of the dynamic boundary condition makes it very straightforward to analyse the motion of the contact point, x = X c (t), (which we have avoided in this formulation by taking > 0) where continuity of both components of the fluid velocity means that U = t · W and γ = 0. Only two terms in (5.6) are therefore non-zero at the contact point, and we obtain consistent with the results of § 4, but valid for all t > 0. This relationship between the acceleration of the contact point and the slope of the free boundary at the point of contact shows that, for non-zero Richardson number, the initial acceleration is infinite when t = 0, with the results of § 4 giving d 2 X c /dt 2 = O(−log t) as t → 0. It also shows that any steady solution must have dY/dX = 0 at the contact point, i.e. the free surface must be tangent to the boundary. We have been unable to construct a non-trivial local solution of this nature, which suggests that steady solutions are ruled out by the behaviour of the solution close to the contact point. In terms of physical variables, abusing our notation for one equation, Finally, when Ri = 0, equation (5.6) shows that dX/dY = 0 at s = 0, i.e. that the free surface is perpendicular to the boundary at the contact point when gravity is negligible. Following Baker & Beale (2004), we define a vortex blob method by modifying (5.5) to be (5.9) We use the lowest-order choice of regularisation function, g(r) = 1 − e −r 2 , as discussed by Baker & Beale (2004), since we are not as concerned with the details of the flow after the curvature singularity has been regularised, and because we use a low-order spatial discretisation. We truncate the surface and discretise using Lagrangian markers and linear elements, with X = X j , Y = Y j , γ = γ j and U = U j at α = j, for j = 1, 2, . . . , N, so that there are 4N unknowns. The integral in (5.9) is evaluated using four point Gaussian quadrature. The parameter α is taken to be gridpoint number, and all spatial derivatives are calculated using central differences. Timestepping is done using the implicit midpoint method, solving the resulting system of nonlinear algebraic equations using quasi-Newton iteration at each timestep. We find that we can reuse the Jacobian, the calculation of which using finite differences is the main numerical bottleneck in the simulation, over multiple timesteps. Note also that the simple, linear representation of the free surface means that only a small portion of the integral in (5.9) has to be recalculated when finding the Jacobian, significantly speeding up the calculation. We chose to use an implicit timestepping method in order to reduce numerical instability and to be able to time step for accuracy rather than stability. In spite of this, for sufficiently large Ri, there is a weak, grid scale instability in γ , which we damp using five point smoothing at each timestep (cf. Longuet-Higgins & Cokelet 1976). The kinematic and dynamic boundary conditions, equations (5.4) and (5.9) provide 3N algebraic equations, and we also specify the required element sizes/marker point separations, s j , using This allows us to control the tangential spacing of the marker points at each timestep, which is possible because we use implicit timestepping and solve for the tangential marker point velocity, U. At each timestep we check whether the curvature of the boundary is adequately resolved, and if not, regrid using quintic splines, so that the curvature of the surface is represented by a cubic spline. We begin our discussion by making a comparison between the asymptotic solution forρ 1, which we discussed in § 3, and the numerical solution of the fully nonlinear free boundary problem.

Comparison of numerical and asymptotic solutions for smallρ
We first consider the solution when Ri = 0 andρ = 0.01. Note that when Ri = 0, there is no initial singularity at the contact point, so we solved the original problem, with no precursor layer ( = 0). Figure 8 shows the position of the contact point as a function of time. The agreement between the asymptotic and numerical solutions is excellent at relatively small times, but has a later discrepancy, which is within the range we would expect given the accuracy of the spatial discretisation of the free boundary in the numerical solution. At times before approximately t = 5, there is good agreement between the numerical and asymptotic solutions in this case, for both the free surface and γ . This can be most clearly seen in figure 9, which shows µ s ≡ γ /s α at various times up to t = 20. The agreement between asymptotic and numerical solutions at earlier times is clear. At later times, µ s steepens, as expected, but in addition, oscillations develop, which, within the framework of a weakly nonlinear theory, indicates the formation of a dispersive shock. It should be possible to develop a weakly nonlinear theory, but we have not attempted this here. The comparison between numerical and asymptotic solutions for µ s is a good test of the numerical method, since the asymptotic solution comes from the analytical solution of the simple, semilinear equation (A 6), whilst the numerical solution for γ comes from solving the evolution equation (5.6) on the evolving free boundary. Figure 10 shows a comparison of the positions of the free surface at later times, when an increasing discrepancy develops between the asymptotic and numerical solutions. This is due not to discretisation error, but to nonlinearity affecting the deepening depression in the free surface. As the slope of the free surface increases, oscillations form at the rear of the depression, which are coupled to the oscillation in µ s that we can see in figure 9. The numerical solution shown was computed using Method 1, and we expect that it will terminate in a curvature singularity at some later time, beyond the reach of our numerical method. We can however see this occurring for slightly largerρ, so we now consider the numerical solution when Ri = 0 and ρ = 0.1. This value is not particularly small, but leads to behaviour qualitatively similar to the solution forρ 1, as shown in figure 11. For this solution, which was computed with δ = 0.05, we can also see the depression in the free surface, followed by steepening and the development of a dispersive shock, but now we can solve at times large enough that overturning and the formation of a Kelvin-Helmholtz roll can clearly be seen. Solutions for non-zero values of Ri (and = 10 −3 ) are shown in figure 12. As Ri increases, the position of the Kelvin-Helmholtz roll moves towards the head of the flow. Note that for small enoughρ, as is the case in figure 12, there is no overturning of the surface at the head of the layer. We should not expect the solutions in figure 12 to agree with the asymptotic solution forρ 1, since this was constructed on the basis that Ri = O(ρ). Figures 13-15 show the effect of varyingρ with Ri fixed at 0.1, 2 and 10 in successive figures. In each case, we again computed until 1500 marker points were needed to resolve the free surface, which was generally just before pinch off. Movies that show the evolution of the surface for these, and other, combinations of Ri andρ are available as Supplementary Material to this paper at https://doi.org/10.1017/jfm. 2019.273, as is the MATLAB code used to implement the adaptive vortex blob method. When Ri = 0.1 (figure 13), and hence the effect of gravity is weak, the general tendency is for Kelvin-Helmholtz instability to occur on the upper surface of the layer, with larger values ofρ leading to a stronger Bernoulli effect, and fluid being projected upwards more strongly. There is no sign of any overturning at the head of the layer when Ri = 0.1. When Ri = 2, and the effect of gravity is stronger ( figure 14) overturning can also occur at the head of the layer, but for sufficiently largeρ, the initial overturning region is transported along the layer, and fluid is again projected upwards. It is worth noting that the solutions shown in figure 14 do bear some resemblance to the two-dimensional large eddy simulation results presented in figure 6 of Carroll et al. (2012), which also suggest that there are large scale bulk motions that are not simply Kelvin-Helmholtz waves. When Ri = 10 (figure 15), and the effect of gravity is even stronger, overturning at the head of the layer is the dominant effect for sufficiently largeρ.

Conclusion
In this paper we have considered the two-dimensional, irrotational flow of a layer of inviscid fluid injected into a different inviscid fluid of lower density by a source on a horizontal boundary moving at constant speed. We were motivated by the work of Carroll et al. (2012), which suggests that this flow provides a very simple model for a powder snow avalanche. With this in mind, after finding that the initial, dambreak dynamics gives rise to complex, multiscale behaviour near the intersection point between the fluid layer and the solid boundary, we chose to regularise this motion by including a thin precursor layer ahead of the original fluid layer. In addition, we chose to deal with the Moore singularities that inevitably arise in an inviscid two fluid flow by using a vortex blob method instead of the more usual surface tension regularisation, which is not appropriate in the context of powder snow avalanches. We also introduced an adaptive vortex blob method. This can either be used to integrate up to the formation of a Moore singularity (Method 1, δ = 0), or, by choosing δ > 0 to be a dimensionless length that characterises small scale processes not modelled by inviscid, irrotational flow, to compute the roll up of the interface on this length scale (Method 2). We found that for small enough dimensionless density difference or Richardson number, the interface rolls up and overturns away from the head of the flow, but that once both the Bernoulli effect and gravity are strong enough, overturning occurs at the head of the layer.
There are some issues that arise from this work, which were mentioned in the paper. Firstly, it would be interesting to derive the weakly nonlinear evolution equation for the free boundary whenρ 1 and t = O(ρ −1 ), which numerical solutions suggest should describe the formation of a dispersive shock. Secondly, the structure of the asymptotic solution when 1, and its relation to the asymptotic solution when = 0 and t 1 is not clear, and is a challenging problem to try to solve using the method of matched asymptotic expansions. Finally, although it is not relevant to avalanche flows, regularisation of the singular solution when = 0 and t 1 by the inclusion of surface tension would be of interest in other two-fluid dam-break problems (see Uddin & Needham (2015) for the single fluid dam-break problem). Other extensions to this work that would be of interest in the context of modelling powder snow avalanches include: unsteady one-dimensional modelling of this flow, perhaps using the ideas described in Khodkar, Nasr-Azadani & Meiburg (2017); the use of a level set method to advance the solution beyond the self-intersection time and overcome the main drawback of the boundary integral method that we have used in this paper; including the effect of an underlying slope; adding feedback between the flow and the speed and strength of the source -this is straightforward to include in the existing numerical method, but a sensible model for this feedback needs to be determined, perhaps based on the work of Louge et al. (2011), Carroll, Louge & Turnbull (2013. In Carroll et al. (2013), a steady, one-dimensional model is proposed, with an influx of fresh particulate material prescribed based on a simple model for the cohesive failure of fresh snowpack under the action of the pressure gradient in the eruption current. The challenge is to use these ideas in the context of an unsteady, two-dimensional flow. Acknowledgements I'd like to thank my colleagues B. Turnbull (whose work funded by Royal Academy of Engineering award RB4525 gave rise to this idea for modelling powder snow avalanches) and M. Scase for bringing this problem to my attention, and also for helpful discussions and comments on a draft of this paper.
We begin by studying the dynamics of the vortex sheet strength, which decouples from the rest of the flow in this limit.
The steady solution of (A 6) is which is shown for various values ofRi in figure 16. We find that Since R 0 ∼ (π − θ) −1 as θ → π, we can write this as which shows that the jump in tangential velocity tends toRi + 1/2 downstream, as x → −∞.
In order to understand the unsteady dynamics of (A 6), note that the solution develops on characteristic curves given by which pass through θ = θ 0 when t = 0, as shown in figure 17. Since R 0 ∼ 1/(π − θ ) and dθ /dt ∼ (π − θ) 2 as θ → π, we find that, in terms of R 0 , the characteristics have dR 0 /dt → 1 as R 0 → ∞, as shown in figure 18. Physically, this means that, away from what we can see in figure 18 is a small neighbourhood of the head of the region occupied by fluid 1 (at x = 1/π, y = 0), disturbances propagate along the interface at unit speed, which is the dimensionless speed of the source. The closer to x = 1/π the characteristic originates, the longer it remains in its neighbourhood before escaping to infinity, with θ 0 ∼ θ e −πt as θ → 0. The solution of (A 6) is then given by  figure 17, redrawn in terms of R 0 (θ). Note that R 0 (θ) is a monotonically increasing function with R 0 (0) = 1/π.
The solution with initial condition µ(θ, 0) = 0 is shown in figure 19 for various values ofRi. It is straightforward to evaluate this accurately for any given values of θ and t by solving (A 10) numerically to find θ 0 and then evaluating (A 11) numerically. The solution can be seen to have two distinct regions: 0 θ < θ i (t), in which the solution is given by the steady form µ = µ 0 (θ), and θ i (t) < θ π, in which the solution grows with t. We can confirm this analytically by considering the large time form of the solution. . The solution, µ(θ, t), when µ(θ, 0) = 0 and t = 0.5, 1, 1.5, 2, 2.5, 3, 3.5 and 4, plotted as a function of θ (a) and R 0 (b), forRi = 0, 1 and 10. The steady state solution is shown as a broken line.
A.2. Solution for t 1 It is clear from (A 10) that when t 1, θ 0 1 and µ ∼ µ 0 , unless θ is close to π. By considering the forms of (A 7) and (A 10) when π − θ = O(t −1 ), we find that the asymptotic solution does indeed consist of two distinct regions, separated by a thin, interior layer.
In summary, as t → ∞, oe/π µ FIGURE 20. The solution, µ(θ, t), when µ(θ, 0) = 0,Ri = 10 and t = 4, along with the asymptotic solution for t 1 in Regions I and III. The position of the interior layer, at θ = θ i , is also indicated. The value of θ i is not very accurate for this low value of t, but for sufficiently large values of t the asymptotic and numerical solutions become indistinguishable.
Note that the solution in Region III shows that the sign of dµ/dθ changes from negative to positive asRi − 1 increases past zero when t 1, which can clearly be seen in figure 19, and indicates a change in the nature of the solution in the far field. The asymptotic solution for t 1 is in excellent agreement with the exact solution, even for low values of t. The solution when t = 4 andRi = 10 is shown in figure 20, along with the asymptotic solution in Regions I and III. Far downstream, as r → ∞, This shows that the asymptotic expansion (A 1) becomes non-uniform when t = O(ρ −1 ) whenρ 1, and we expect that the solution becomes nonlinear on this time scale. We will discuss this further in § 5.
A.3. The normal velocity at the free surface, n · ∇φ 1 = n · ∇φ 2 In order to determine the normal velocity at the free surface, which controls its motion, we use the vortex sheet representation (see, for example, Baker et al. 1982), where t(s) ≡ (X s (s), Y s (s)) is the unit vector tangent to the boundary and µ s ≡ ∂µ/∂s is the jump in the tangential velocity at the free surface, also known as the normalised vortex sheet strength. We have extended the free surface by reflection in the x-axis, with µ s (−θ , t) = µ s (θ, t), X(−s) = X(s) and Y(−s) = −Y(s) by symmetry. The arclength along the fixed boundary can be calculated from dθ ds = π sin 2 θ sin 2 θ − 2θ sin θ cos θ + θ 2 .

(A 19)
Equation (A 18) allows us to determine the normal velocity from the vortex sheet strength, µ s , which we discussed above. Note that it is possible to use conformal mapping to find the individual velocity potentials,φ 1 andφ 2 , as integrals of n · ∇φ i , but we will not present any details here as it is possible to determine the motion of the free surface from the normal velocity alone. Figure 21 shows how the vortex sheet strength varies for the same values of time, t, used in figure 19. We can see that µ s approaches the steady state solution in Region I as t increases, following the asymptotic solution described in § A.2 even at moderate times. In asymptotic Region III, where R > t − (1/π)(log t − 1) + O(t −1 ), In order to determine the behaviour of the free surface, we use the known functional forms of φ 0 (r, θ ) and R 0 (θ) to show that the radial perturbation of the free surface satisfies In order to solve the simple, linear equation (A 22) numerically, we discretise in θ and approximate ∂R/∂θ using first-order upwind finite differences. The resulting system of ordinary differential equations can then be solved using ode45 in MATLAB, coupled with the simple initial value problem for θ 0 (t; θ ), This follows from integrating back along a characteristic at fixed θ, and allows θ 0 , and hence µ s , to be calculated efficiently as part of the solution. The principal value integral in (A 18) can be evaluated numerically, working in terms of θ, subtracting out the singular behaviour in the neighbourhood of s = s (θ = θ ) and using four point Gaussian quadrature. This linear relationship means that the discretised version of n · ∇φ i can be found from µ s (θ; θ 0 (t)) by multiplication by a constant matrix, which can be precalculated. Some typical solutions are shown in figures 2 and 3. The asymptotic solution for t 1 can be constructed in terms of an inner and an outer asymptotic region.
B.1. Outer solution: At leading order, we obtain the boundary value problem Note thatX(y) decouples from the outer problem, and can be determined from (B 5) once the normal velocity at the free surface is known.
B.1.1. Local solution as (x − 1/π) 2 + y 2 → 0 Close to the point x = 1/π, y = 0, where the undisturbed free surface meets the boundary (a stagnation point of the initial potential, φ 0 (x, y)), the incompatibility in the boundary conditions (B 3) and (B 4) leads to a logarithmic singularity, as discussed at the end of § 3. The local solution is {r log r cos θ − θr sin θ + kr cos θ} as r → 0, (B 7) where (r, θ ) is here a set of polar coordinates centred on (1/π, 0), and p 0 and k are parameters that are determined globally. Although p 0 ≡ p 0 (ρ,Ri) is unimportant, k ≡ k(ρ, Ri) affects the solution in the inner region. By using (B 7) in (B 5), we find thatX It is this logarithmic singularity in the position of the free surface in the outer solution that leads to an inner region where the singularity is resolved. We will also use the streamfunctions,ψ 1 andψ 2 , that correspond to the velocity potentialsφ 1 andφ 2 , which havē {r log r sin θ − (π − θ )r cos θ + kr sin θ} as r → 0, (B 9) {r log r sin θ + θ r cos θ + kr sin θ} as r → 0.
(B 10) We will use capital letters with subscripts to denote the potentials and streamfunctions on the boundary as a function of arclength, s, and note that (B 12)

B.1.2. Solution using Plemelj formulae
We define the complex potential, Φ(z) = φ i + iψ i in D i , with z = x + iy, which can be written in terms of the potential difference on the boundary, ∂D, which we extend into the lower half-plane by symmetry, as where Z 0 (s) ≡ X 0 (s) + iY 0 (s). On the boundary, the Plemelj formulae (see, for example, Ablowitz & Fokas 1997) show that We can now use (B 4) to eliminate Φ 2 , and obtain a single Fredholm integral equation of the second kind, Note that the kernel is regular as s → s, so this is not a principal value integral. When ρ = 0, this immediately gives Φ 1 (s), and forρ > 0 is in a form that can be solved iteratively, using four point Gaussian quadrature to evaluate the integral. In order to determine the displacement of the boundary,X, note that the streamfunction on the boundary is After regularising this integral and subtracting out the behaviour of Φ 1 (s ) close to s = 0, which leads to the logarithmic behaviour of Ψ as s → 0, equation (B 12), we can evaluate Ψ using four point Gaussian quadrature. We can calculate ∂Ψ /∂s = n · ∇φ i , and hence the displacement of the boundary,X, using (B 5). Figure 22 shows the solution for various values of Ri andρ, although we have not plotted a point at y = 0. As expected, there is a logarithmic singularity inX(y) as y → 0, but it is too weak to be seen clearly in figure 22. The head of fluid 1 swells and is displaced forwards close to the contact point, reproducing the behaviour that we saw forρ 1 andRi = O(1) in figure 2. The rate of displacement and size of the swelling increases withρ and Ri, but is always qualitatively similar.
B.2. Inner solution: φ i = O(t 3 (−log t) 2 ), x − 1/π, y = O(t 2 (−log t)) On comparing the size of the terms neglected in the outer region to those retained as s → 0 (equivalently, y → 0 on the boundary), we find that appropriate scaled variables are (see Uddin & Needham (2015) for more details on the very similar, single fluid version of this dam-break problem) x = 1 π + Rit 2 (−log t)x, y = Rit 2 (−log t)ŷ, X = 1 π + Rit 2 (−log t)X, (B 17a−c) and matching with the outer solution requireŝ This suggests that we should expand the unknowns aŝ On substituting these expansions into (B 19) and (B 20) we find that, although we obtain a fully nonlinear free boundary problem at leading order, there is a simple, linear solution that reproduces the matching conditions up to a constant. A simple linear solution is also available at O(log(−log t)/(−log t)), and we find that The constants K i and L i are not fully determined at this order.
At O((−log t) −1 ) we obtain a non-trivial boundary value problem set in the quarter planesŷ > 0 andx <X 0 forφ 12 andx >X 0 forφ 22 . After definingφ i ≡φ i2 ,X ≡X 2 , x =x −X 0 andỹ =ŷ, we arrive at∇ Note that, at this order in (−log t) −1 , the inner region does not see the leading-order flow, φ 0 (x, y), in these matching conditions, since φ 0 is quadratic in r as r → 0 close to the stagnation point. This is why this inner problem is equivalent to that for a two-fluid dam-break problem, with initially stationary fluids.
A key feature of this boundary value problem is the term involvingρX on the righthand side of (B 30). This is the reason that the limitρ → 0 is singular. It is this term, which arises from the inertial terms in the Bernoulli equation in this accelerating inner region, that allows a bounded solution to exist. Simply differentiating (B 30) with respect toỹ, taking the limitỹ → 0 and using (B 31) shows that ∂X ∂ỹ (0) = − π 1 + 1 2ρ 2ρ .

(B 35)
In terms of the original variables, this is ∂X ∂y (0) ∼ − π 1 + 1 2ρ 2ρ(−log t) as t → 0, (B 36) which reproduces the result of King & Needham (1994) and Uddin & Needham (2015) in the single fluid limit,ρ → ∞. Note that it is also possible to recover (3.28)-(3.31) with σ = 1 in Uddin & Needham (2015) in this limit after a rescaling. In the Boussinesq limit,ρ → 0 with Ri = O(1), the slope at the contact point and, as we shall see, its position, tend to infinity, and hence this limit is a singular perturbation.
produced what appears to be the same oscillatory solution (figure 4). In Needham, Chamberlain & Billingham (2008, appendix C), similar oscillations were shown to exist close to an inclined accelerating plate. A local analysis of (B 43)-(B 46) above leads to an equation that controls the behaviour of the oscillations which is identical to (C.9) in Needham et al. (2008, appendix C) with γ = 2, consistent with the asymptotic expansion (B 56). Note that, although ξ → 0 and ξ Y → 0 as Y → 0, the curvature of the free boundary is singular, with ξ YY = O(Y −3/2 ) as Y → 0. In terms of the original variables, we havẽ X = 1 2π 1 + 1 2ρ 2 π 1/2 Γ (B 58) The logarithmic singularity in this expression asρ → 0 indicates that this is a singular perturbation.