Axisymmetric plumes due to fluid injection through a small source in a wet porous medium

A small spherical source discharges a fluid into a porous medium that is already fully saturated with another fluid. The injected fluid has higher density than the ambient fluid, and so it forms a plume that moves downward under the effects of gravity. We present a simple asymptotic analysis assuming the two fluids do not mix that gives the width of the plume far from the source as a function of the injected volume flux. A spectral method is then developed for solving the full nonlinear problem in Boussinesq theory. Accurate numerical solutions are presented, which show in detail the evolution of the plume of heavier injected fluid over time. Close agreement with the asymptotic plume shape far from the source is demonstrated at later times.


Introduction
The study of the flow of fluid through rock or soil, and the path it takes, is relevant to many practical scenarios. An example of such a flow arises in the contamination of groundwater by an industrial leak [1]. Further applications arise in in situ mineral leaching, in which a lixiviant is injected into rock to dissolve the ore and recover it in solution [2,3], and the extraction of fluid from groundwater aquifers or oil reservoirs [4,5].
Flow in fully saturated porous media is generally modelled using a well-known system of linear differential equations. When phreatic surfaces or interfaces are present, however, the location of these free boundaries is usually unknown. This makes the mathematical problem highly nonlinear, and thus, difficult to solve, since both the B Catherine A. Browne catherine.browne@utas.edu.au B Lawrence K. Forbes larry.forbes@utas.edu.au 1 School of Natural Sciences, University of Tasmania, Hobart, TAS 7005, Australia fluid velocity and the shape of the fluid region itself must be determined as part of the solution process. Closed-form solutions to these problems, in terms of familiar analytical functions, are, therefore, not usually possible, although some classical results are known for problems in which a high degree of symmetry exists and certain simplifying assumptions are possible (see the texts by Bear [6,7] and Strack [8]). Otherwise, numerical methods must be used instead.
Many of the studies of plumes in porous media involve an assumption that the flow can be regarded as two dimensional. In addition, by considering the flow that is established after a long time, only steady-state behaviour need be considered, and this simplifies the problem description substantially. This has allowed a number of important insights into plume structure to be obtained. Wooding [9] carried out a twodimensional steady asymptotic analysis in which the outer fluid could be entrained into the moving plume, which as a result becomes wider as it develops. This analysis is valid when the vertical speed in the plume is almost equal to that in the ambient fluid, and predicts that the width at the head of the entraining plume increases as the 2/3 power of its distance from the source. In an experimental investigation, Oostrom et al. [10] found that leachate plumes may be stable or unstable, depending on the density ratio of the two fluids and the volume flow rate at the source, and conditions for unstable plumes to form fingers were studied numerically and experimentally by Cremer and Graf [11]. Hewitt et al. [12] considered a steady two-dimensional model for a buoyancy-driven plume that encounters an internal layer of different permeability, and carried out a detailed asymptotic analysis of conditions under which the plume spreads. A very detailed numerical analysis of two-dimensional plumes has been undertaken by Slim [13], in a study that now incorporates time-dependent effects. She finds six different time regimes for plume evolution, including an early stage during which the plume widens and an intermediate stage during which dense fluid is stripped from the plume interface as it accelerates downwards.
In a recent paper, Browne and Forbes [14] analysed the formation and development, through time, of a vertical plume formed by the injection of a more-dense fluid into a wet porous medium. They considered two-dimensional geometry only, in order to simplify the mathematics. Thus, the flow occurred in vertical planes, with the heavier fluid injected through a horizontal line source orthogonal to these flow planes. Horizontal drilling wells are, in fact, used in some mining applications [15,16]. Browne and Forbes [14] solved their planar flow problem using a spectral approach, in which the fluid velocity and density were represented using Fourier series with time-dependent coefficients. This gave highly accurate results, and the evolution of the fluid region from an initial cylinder into a fully developed vertical plume was demonstrated.
The purpose of this present paper is to investigate the analogous problem, in which the flow geometry is now three-dimensional and axisymmetric, with the denser fluid injected through a small spherical source. This would represent a situation as occurs in mineral leaching, for example, in which the leaching liquor is introduced through a narrow vertical pipe drilled into the rock. The equations describing this flow are presented in Sect. 2. To solve this nonlinear problem numerically, we use a spectral method similar in some ways to that in Browne and Forbes [14], but here we develop new basis functions, as discussed in Sect. 3 that allow us to consider a finite-sized source region. This is possibly more realistic and closer to industrial reality than approximating the injection region as an idealized point source, and it avoids the singularity that a point source represents. Nevertheless, at some distance from the source, outflow from a point source or a finite source would be almost indistinguishable from one another [17], as also observed experimentally by Bietenholz et al. [18]. We also present an asymptotic closed-form solution in Sect. 4 that estimates the plume shape far from the injection zone, and we compare this approximate solution with the computations of the fully nonlinear problem in Sect. 5. We consider the timedependent behaviour of this three-dimensional flow, and a key aim is to determine the way in which the plume evolves in time, towards a steady-state configuration. A discussion in Sect. 6 concludes this paper.

Governing model
We wish to model the trajectory of some fluid as it flows through a porous medium of uniform porosity n and permeability k. The fluid is injected from a point source at the origin of a Cartesian coordinate system where the z-axis points vertically and the horizontal ground is located parallel to the x y-plane high enough above the point of injection such that the fluid remains underground. The acceleration of gravity, g, is directed downward. We assume that the fluid has density ρ 1 , and that the porous rock is already saturated with a fluid of a lower density, ρ 2 ; however, both fluids have the same constant viscosity, μ.
Fluid is flowing radially in all directions from the injection point, and given this geometry, it is intuitive to use a spherical coordinate system, which will be done for the rest of this paper. We assume that the flow of fluid is axisymmetric such that there is no dependence on the azimuthal angle θ , but only the radial distance r and the polar angle φ measured from the z-axis.
Initially, we assume there is a sharp interface between the region the heavier fluid saturates and the rest of the porous medium. As fluid is injected into the system, the two fluids will mix and it is no longer appropriate to assume a sharp boundary between two distinct fluids. We instead assume that the fluid density ρ(r , φ, t) varies continuously across a narrow boundary of finite width from ρ 1 (the density of fluid injected at the source) to ρ 2 (the density of the fluid that saturates the ground far away from the source). Similarly to a number of other papers [14,[19][20][21], we assume that the higher density of the introduced fluid is due to the addition of some solute with concentration C(r , φ, t) such that the density varies linearly with the concentration and C → C 0 as ρ → ρ 1 ; that is The concentration of the solute satisfies the transport equation: where D is a constant diffusion coefficient for the solute and n is the porosity of the medium, a dimensionless constant that gives a measure of the void space in the porous rock that the fluid can flow through. From the pore pressure of the fluid, p, we can find the seepage velocity q using Darcy's law in the form [22], Here μ is the dynamic viscosity of the fluid and k is the permeability of the porous medium. The unit vectors e r and e φ point in the directions of increasing radial r and angular φ coordinates, respectively. We use a Boussinesq approximation in which the fluid velocity also satisfies the continuity equation: representing an incompressible fluid. We assume that initially, the higher density fluid saturates a spherical region of radius a around the point source and that at time t = 0 the point source is suddenly turned on, thereafter, injecting fluid at the volumetric flow rate Q. The velocity of the fluid near the source can be written as follows: At this point, it is convenient to make use of dimensionless variables, where all lengths have been scaled by the initial radius a of the injected fluid. Similarly, all times are scaled by nμa/kρ 2 g, pressures by ρ 2 ga, and the concentration of the solute is scaled with reference to its value C 0 at the source. Using this scaling, our governing Eqs. (1)-(3) become and The non-dimensional velocity of the fluid near the source is given by In the equations above, three dimensionless parameters have appeared. There is a density ratio, a Froude number, and a dimensionless diffusion coefficient for the solute, At the point of injection, the scaled concentration of the solute is subject to the boundary condition that C → 1 as r → 0.
Far from the source, Given that the velocity of the fluid satisfies the continuity Eq. (4), we know that a stream function (r , φ, t) exists such that the fluid seepage velocity can be written as follows: Hence, the condition (9) for the velocity near the source can be transformed into a condition on the stream function; that is, Using the state Eq. (6) and the stream function form for the velocity as in (15), by taking the curl of Darcy's law (8), we can show that

Numerical method
In this section, we describe the numerical method we use to solve (17) coupled with the transport Eq. (7), along with the conditions (16) and (13) at the source. It follows the method given in Browne and Forbes [14] which solved a similar system of equations, but for planar flow involving fluid injected along a line source. It is necessary to create an artificial boundary at r = R ∞ that is chosen to be large enough that it does not interact with the top of the plume; naturally it is inevitable that the bottom of the plume will always reach the computational boundary given enough time. At this outer boundary, we follow (14) and introduce the condition that meaning that there is no solute in the fluid far away from the source. We have also introduced a boundary at r = α around the source to exclude the origin from the computational volume, which consequently is now a spherical shell as shown in Fig. 1. This was done to avoid problems caused by the singularity at the origin, but in any event, the finite fluid speed on the injection sphere r = α more closely mimics physical reality. We replace the boundary conditions (13) and (16) at the source with conditions at this inner boundary such that and Note that the computational boundaries are chosen such that The dimensionless concentration of the solute is represented in the form: where p n+ 1 2 (κ m,n r ;α) is the Bessel function cross-product and κ m,n is the m th positive root of Here, J n+1/2 (z) and Y n+1/2 (z) are the Bessel functions of half-integer order of the first and second kinds, respectively, and we have used the notation p n+ 1 2 (κ m,n r ;α) for the Bessel function cross-product, so as to be consistent with Abramowitz and Stegun [23, p. 361]. We can numerically find these roots in (24) using the method as described by Horsley [24,25]. Not to be confused with p n+ 1 2 (κ m,n r ;α), the terms P n (cos φ) are the Legendre polynomials of order n, as defined in Abramowitz and Stegun [23]. We note that the basis functions go to zero on both boundaries, so the boundary conditions (18) and (19) for the concentration of the solute are satisfied by the first two terms in (22).
The stream function is also represented spectrally: the form ensures that (17) is satisfied when the B m,n (t) are chosen appropriately. In these series representations for C(r , φ, t) and (r , φ, t), the coefficients A m,0 (t), A m,n (t), and B m,n (t) are unknown functions of time t and the expressions become exact as the numbers M and N of Fourier modes approaches infinity. The first term in (25) is chosen to satisfy the condition (20) that models the fluid being injected into the system and the second term is there to represent the downward acceleration of gravity. As they satisfy they do not have equivalent terms in the series for the concentration (22). For the B m,n (t) to satisfy equation (17), we must find the relation to the A m,n (t) using Fourier analyses. We substitute the spectral forms for the concentration (22) and (25) into (17), and after some algebra, we find that (κ m,n r )P n (cos φ) sin φ We then multiply (27) by functions r 3/2 p l+ 1 2 (κ k,l r )P l (cos φ) sin 2 φ and integrate over the volume α < r < R ∞ , 0 < φ < π to get where and Note that where we have used the same notation for the derivative of the cross-product as in Abramowitz and Stegun [23, p. 361].
To solve numerically for C(r , φ, t) (and consequently (r , φ, t) as well), we wish to find the coefficient functions A m,0 (t) and A m,n (t) so as to satisfy the transport Eq. (7) for the concentration of the solute. The coefficient functions B m,n (t) will follow using (28). We use Fourier analysis to find the following system of nonlinear ordinary differential equations for the coefficients: and We also need to transform the initial conditions for the concentration, C(r , φ, 0) into initial conditions for the coefficients, A k,0 (0) and A k,l (0). Recall that initially we assumed that the porous medium is already saturated with some fluid without the presence of the solute (corresponding to a concentration of the solute equal to zero) and that there is a spherical region with radius r = r init around the source in which a solute is dissolved in the fluid such that the concentration of the solute is equal to one. Therefore, where r init is another parameter to be chosen. Using similar Fourier analysis on our series representation (22), we can show that the initial coefficients are and We make use of Lanczos smoothing [26] to lessen the impact of Gibbs' phenomenon [27] occurring in our spectral solution at the jump discontinuity at r = r init . This has been done by multiplying each Fourier coefficient A k,0 (0) by sin( k)/ k where is a small parameter. We typically chose this Lanczos parameter to have the value = 0.1.
The system of Eqs. (32) and (33) is then integrated numerically forward in time using the MATLAB routine ode45. To do this, we have evaluated the integrals in (32) and (33)

Asymptotic solution
As a comparison to our numerical results, we also look for a large-time steady-state solution by considering the scenario in which there exists a sharp interface r = R(φ) between the two fluids. The densities are constant either side of this interface and the fluids cannot cross this sharp boundary. As such, and turning our focus to the denser fluid being injected, the fluid is subject to the kinematic condition: where u and w are the velocity components such that q = ue r + we φ . As we are looking for a steady-state solution, we want to find an interface R(φ, t) ≡ R(φ) that does not depend on time. We assume the approximate asymptotic form: where the first term represents the outflow of fluid from the point source as in (16) and the second term represents the downward acceleration of gravity. It implies that u = F 4πr 2 − cos φ and w = sin φ.
Thus, far from the source at the bottom of the plume (i.e. as r → ∞ and φ → π ), the velocity of the descending plume approaches q = 1e r in these dimensionless variables. Substituting this solution for into the kinematic boundary condition and letting ∂ R/∂t = 0, we find that This is a first-order differential equation for R(φ) which we can solve to find where γ is a constant that is yet to be determined. We choose γ by requiring R(φ) to be finite at the top of the plume where φ = 0 and as such, choose γ = F/2π . Therefore, as time approaches infinity, we find that the asymptotic interface becomes Fig. 2 Plots of the initial concentration profile as defined in (34) for r init = 1 on the left, and r init = 0.45 √ F/π for F = 5 on the right. The profile has been reconstructed from the Fourier representation (22), and for both plots, we have used M = N = 81 Fourier modes over the domain 0.3 < r < 10 with a Lanczos smoothing parameter of = 0.1. In both plots, evidence of Gibbs' phenomenon [27] is visible as small wavelets near the jump discontinuity In Sect. 5 below, we compare this asymptotic solution to the results of our numerical solution. We also use it to guide our choice of r init , the radius of the initial spherical region inside which C = 1 at t = 0. From (42), we can show that far from the source (i.e. as φ → π and r → ∞) the radius of the plume approaches √ F/π. For the results shown in this paper, we have set r init as the smaller of one and 0.45 √ F/π, both of which are illustrated in Fig. 2. We have done this to balance between avoiding an initial condition that is larger than the asymptotic plume shape and also avoiding an initial condition that fills too much of the computational domain.

Results
In this section, we present a number of numerical results for a variety of different parameter values. Recall that our three-dimensionless parameters in (10)-(12) are a density ratio D, a Froude number F, and a dimensionless diffusion coefficient β. We have used D = 1.05 for all results shown in this paper, which is well within the range for the Boussinesq approximation to hold true. We will present results for F = 5, 50 and 100 and β = 0.01 and 0.0001, and we have chosen these to be consistent with the dimensional parameters as follows: •  The above dimensional parameters give rise to a large range of values for F and β, and given the constraints of time and computing power, we have only found numerical solutions for a subset of sensible values. However, the results shown in this paper are consistent with the physically realistic dimensional parameters as described above.
In order to establish the convergence of the numerical solution, it is important to demonstrate that the numerical results are independent of the numbers M and N of Fourier modes in our series representations (22), (25). In Fig. 3, we have compared colour maps of the concentration where we have used the same parameter values in all, except for the number of modes and mesh points used. The picture on the left in Fig. 3 was computed using M = N = 31 Fourier modes with 155 × 155 grid points over the computational domain. Although evidence of Gibbs' phenomenon can be seen over the whole domain, the shape of the descended plume can still be sufficiently observed. In comparison, the middle diagram of Fig. 3 was computed using M = N = 61 with 305 × 305 grid points over the computational domain and shows a more converged numerical result. The shape of the plume is now much sharper and although still present, oscillations due to Gibbs' phenomenon are greatly reduced. Lastly, the rightmost picture in Fig. 3 was computed using M = N = 81 Fourier modes with 405×405 grid point and shows a very converged numerical result with oscillations associated with Gibbs' phenomenon even further reduced. Nevertheless, its plume shape is nearly Fig. 4 Colour maps of the concentration of the solute C for the Boussinesq solution at the dimensionless time t = 6. The solutions were computed using Froude number F = 50, density ratio D = 1.05, and dimensionless diffusion coefficient β = 0.01 using M = N = 81 Fourier modes with 405 × 405 grid points. The solution on the left was solved over the computational domain 0.3 < r < 10, whereas the solution on the right has a larger inner boundary, specifically, its domain is 0.9 < r < 10 indistinguishable from the results with 61 coefficients, showing that convergence has occurred with 81 coefficients. This is consistent with the convergence analysis by Allwright et al. [32]. We note that for Fig. 3 and subsequent figures in this paper, we have shown only a portion of the computational domain such that x ∈ [−6, 6] and y ∈ [− 10,5]. This is the region in which the higher density fluid descends and is therefore the region of most interest. Figure 4 shows a comparison between two solutions computed using Froude number F = 50, diffusion coefficient β = 0.01, and density ratio D = 1.05 at time t = 6 with different inner computational boundaries. The solution on the left was computed over the computational domain from r = α = 0.3 to r = R ∞ = 10 whereas the solution on the right was computed from r = α = 0.9 to the same r = R ∞ = 10. The larger inner boundary means that there is less of a difference between the inner and outer boundaries which makes the solution more numerically stable in respect to the calculation of the roots of the basis functions (23); however, in practice, there is not much difference between the two solutions. Importantly, the shape of the plume is unchanged.
Solutions using Froude number F = 5 are illustrated in Fig. 5 at the three dimensionless times t = 2.75, 5.  the thick black-dashed line. This is independent of β and so is the same for both top and bottom solutions. It can be seen that the steady-state solution is in close agreement with the numerical solutions for both values of β. In both t = 5.5 plots, some disturbance is visible above the top of the plume; by t = 8.25, it can be seen that this disturbance has moved inside of the top part of plume. Nevertheless, it does not appear to have impacted the overall shape of the plume. In addition, at this time, the plume has almost reached the outside edge of the computational domain, thus, running our solution much further past this time has no practical motivation. Initially, at time t = 0, the heavier fluid occupies an initial sphere of radius 0.45 √ 5/π as shown in the plot on the right in Fig. 2, noting that the heavier fluid is in regions where C = 1.
At the earliest time t = 2.75 shown in Fig. 5, the plume widens slightly near its top in both pictures, to form an approximately conical shape. This is caused by the injection of fluid through the finite-radius hole shown, but may also be due to some entrainment of the surrounding stationary fluid in the wetted rock; this is suggested by the paper of Wooding [9] even though the conditions required for that approximation to be valid do not hold here. At later times, however, the plume develops a constant diameter, consistent with the analysis of Slim [13]. Figure 5 shows that as time progresses, the plume then descends under the effects of gravity and fluid injection at the source (9), and given how closely the numerical plume shape follows the large-time asymptotic  Figs. 6 and 7, initially, at time t = 0, the heavier fluid occupies an initial sphere of unit radius as shown in the plot on the left in Fig. 2. We note that the higher the Froude number, the sooner the plume reaches the bottom of the computational domain. This is unsurprising as a larger Froude number (11) corresponds to a larger flow rate Q of the injected fluid (assuming all other parameters are kept constant) meaning that the fluid is entering the system and falling at a faster rate. In addition, the width of the plume widens as the Froude number increases, as suggested by the asymptotic solution (42).

Conclusions
In this paper, we have studied the injection of fluid into a porous medium in threedimensional geometry. We have presented an asymptotic solution, valid for large times, and numerical results using a spectral method. The agreement between the two is very good for later times for all parameter sets we have considered. As the injection rate is increased (representing increasing Froude numbers), the plume that forms becomes wider, and also falls at a faster rate. For plumes with a larger dimensionless diffusion coefficient, the narrow interfacial zone between the heavier and lighter fluids widens and is less sharp. For future work, we may be able to generalize the spectral formulation such that it is able to cope with the presence of multiple singular injection points and extraction points as well.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.