Extraction of density-layered fluid from a porous medium

We consider axisymmetric flow towards a point sink from a stratified fluid in a vertically confined aquifer. We present two approaches to solve the equations of flow for the linear density gradient case. Firstly, a series method results in an eigenfunction expansion in Whittaker functions. The second method is a finite difference method. Comparison of the two methods verifies the finite difference method is accurate, so that more complicated nonlinear, density stratification can be considered. Interesting results for the case where the density stratification changes from linear to almost two-layer are presented, showing that in the nonlinear stratification case, there are certain values of flow rate for which a steady solution does not occur. A spectral method is then implemented to consider cases in which there is a stagnant region beneath a sharp interface between two layers of different, but constant, density. In this situation, flows also exist only for flow rates beneath a critical flux value, consistent with the results for the continuous density stratification.

interfaces or free boundaries using hodograph and boundary element techniques. In all cases, it was found that there was a limit on the steady flow rate beyond which single-layered flows could not be found. This limiting maximal flow is known as coning [14], in which the fluid surface (or interface in the case of a two-layer fluid) is drawn up sharply and enters the sink directly. List [15] and Yih [16] showed that there is a distinct layering of the flow for linear density stratification. At flow rates higher than the critical coning value, solutions include flows with multiple flowing layers, such as those described in [5,[17][18][19]. These solutions are analogous to similar supercritical flows in surface water hydrodynamics [20,21].
Here we examine the case of withdrawal through a point sink in a fluid that has either a continuous density stratification or a distinct layering, and consider the limit as one approaches the other. An eigenfunction expansion [16] is used to solve for the flow if the stratification is linear, and the results are used to verify a finite difference method that is then used to consider nonlinear density strata. After that a spectral method is used to examine the flow into a point sink from a single layer of constant density, separated from a second layer by a sharp density interface. A brief investigation of this problem with the point sink situated on the top of the region was considered in Jose et al. [22] and this work is extended here to consider all elevations of the sink.

Problem formulation
Consider the axisymmetric flow of a density stratified fluid into a point sink within a vertically confined aquifer. In cylindrical coordinates, the sink is situated atr = 0 and at some height,ẑ = H s , within a vertically confined region with height d (see Fig. 1) and coordinates (r ,ẑ). The flow is assumed to be radially symmetric towards the sink, and the fluid to have a densityρ =ρ(r ,ẑ).
Following the work of Yih [16], consider the cylindrical coordinates (r ,φ,ẑ) withẑ increasing vertically upward. In axisymmetric flow, the flow is a function of (r ,ẑ) only, and the pseudo-velocity (u , w ) is defined by where μ 0 is a constant reference viscosity, andû andŵ are the velocity components in ther andẑ directions, respectively. In many situations, with a single fluid, the viscosity μ will be approximately constant at μ 0 , so that the actual velocity and pseudo-velocity will be the same. It is only if the fluid has different properties that this transformation is required. In terms of the pseudo-velocity, the equations of motion invoking Darcy's law become μ 0 Neglecting diffusive effects, the viscosity μ does not change along a streamline in steady flow. The pseudovelocities u and w satisfy the continuity equation and so we define a stream function ψ which is defined by u = − 1 r ∂ψ ∂ẑ and w = 1 r ∂ψ ∂r .
(3) Fig. 1 Sketch defining the variables for flow into a point sink located at (0, H s ). Arrows indicate the direction of flow andρ is the variable density gradient. The region has total height d. The radial coordinate isr and the vertical coordinate isẑ After some work [16], we reach the equation If the fluid is incompressible, then in steady flowρ is a function of ψ only. We now solve (4) together with appropriate boundary conditions where Q is the discharge into the sink. The discontinuity of ψ at r = 0 and z = H s represents the sink. We define the upstream (r → ∞) density asρ → f (ẑ), where f (ẑ) defines the density as a function of height far upstream. To reduce the number of parameters in the problem, we use the dimensionless quantities in which case (4) becomes with the boundary conditions where h s = H s d . If h s = 1, the sink is located on the top surface.

Linear density profile
A long way upstream (ξ → ∞) we assume a density profile of the form ρ = f (ζ ), so that dρ dζ = f (ζ ) is the density gradient. In the case of a linear gradient to be considered in this section, ρ = f (ζ ) = 1 − βζ is the appropriate density function and Ψ → ζ as the flow velocity approaches zero. In incompressible flow, the density along each streamline remains the same, so that it is true everywhere that The nondimensional form of Eq. (4) with the linear density profile is where is proportional to the density gradient or the inverse of the flux into the sink Q. A large value of the nondimensional parameter λ means either a strong density stratification or a low rate of flow into the outlet. The quantity λ is the same as that defined by Yih [16] and we will use it here for consistency. However, we note that if we define a new quantity γ − 1 = ρ 1 /ρ 0 − 1 where ρ 0 is the density at the top of the region and ρ 1 is at the bottom, then β = γ − 1 and we can define a slightly more general quantity where g = g(γ − 1) is a reduced gravity. The quantities λ and G relate the bouyancy of the fluid to the suction at the sink. Thus, the problem is to solve (8) with the boundary conditions (7) along with Ψ → ζ as ξ → ∞. The quantity G will allow us to compare with flows later in this work. In the case of a linear density profile G = λ.

Series solutions
Following Yih [16], and taking h s = 1, i.e. the sink on the top surface, and in accordance with boundary conditions (7), we seek a solution as an eigenfunction expansion in the form The differential equation satisfied by g n (ξ ) is where g n (0) = 0, and g n → 0 as n → ∞. The differential equation (12) has a regular singular point at ξ = 0 and so we use the Frobenius method in powers of ξ [23]. Following Yih [16], we let g n = e σ/2 h(σ ) where σ = λ 2 ξ 2 2 , and then where k n = −n 2 π 2 2λ 2 , n = 1, 2, 3 . . . . Equation (13) has solutions in the form of a Whittaker function [24], and the appropriate solution is At ξ = 0, g n (0) = 1 from (14), and coefficients A n are chosen to satisfy (7), so that the stream function Equally spaced streamlines for λ = 1 and λ = 8 are shown in Fig. 2. For λ = 1, the flow is stronger nearer to the sink due to the higher flow rate, and so the streamlines 'fill' more of the channel. The flow coming from the bottom part of the region is relatively slow when λ = 8 and so the streamlines fill less of the depth. Thus, near the bottom, a long way from the sink, the horizontal flow is much slower. This effect is clearly seen in the velocity profiles shown in Fig. 3. For the case λ = 8, there is a region near the bottom of the aquifer that is almost stagnant. Note that for a linear density profile upstream, the streamlines are also contours of density, and so the local stratification can be seen from the location of the streamlines. Therefore, as λ increases, the density stratification has the effect of restricting the flow so that it is much greater close to the level of the sink. In what follows we will investigate this further for cases in which the point sink is no longer on the top of the region and also for the situation in which the stratification is not linear but takes on a more layered structure as is found in oil reservoirs or coastal aquifers.

Finite difference solution
To compute solutions for flows in which the density stratification is not linear or the sink is located in a different location, we approach the problem using the finite difference method. To proceed, we will compare the solution with that computed by Yih [16] for the case of a linear density variation with the sink on top of the region, presented in Section 3.1, and then consider cases in which the sink is in an arbitrary location, and the density profile is not linear. In particular we are interested in cases in which the stratification approaches two layers of different density separated by a thin interface.
If we make a guess for Ψ (ξ, ζ ), at the grid points (ξ i , ζ j ), i = 2, 3, . . . , N R − 1 and j = 2, 3, . . . , N Z − 1, separated by Δξ and Δζ , respectively, then the errors, E i, j , at the corresponding points for the finite difference form of (8) are The values of Ψ i,1 , Ψ i,N Z and Ψ 1, j are given by the boundary conditions (7). This gives

Changing sink elevation
Henceforth we will assume the sink location to be more general at z = h s as indicated by the boundary conditions (7). Since the sink is now below the top surface, the nondimensional quantity G = λ is changed very slightly because the flux into the sink is now doubled (since the full sink is exposed to the fluid, rather than just the lower half). Therefore, G is now given by to reflect this doubling of the inward flux. The finite difference method was modified to solve the problem with a linear density gradient for arbitrary sink height and the results are shown in Figs. 4, 5, and 6. In Fig. 4, we can see that when the withdrawal rate is larger the streamlines pull closer towards the sink, reflecting the velocity profiles seen in Fig. 5, in which the higher withdrawal rate gives a more even spread over the vertical region. In the case of linear stratification, the streamlines are also contours of density, and so the contours are also evenly spaced in density. When λ = 8, for low flow, the velocity   profile is more confined to a narrow layer in the middle of the region. Note that the quantity λ is proportional to the stratification value β. Thus large λ also corresponds to stronger stratification. This means that for a given flow rate, a stronger density gradient restricts the vertical movement of the fluid leading to the narrow flowing "layer". This layering was identified and studied by List [15] and Yih [16]. In cases where the sink is not at the middle height, i.e. h s = 0.5, as shown in Fig. 6, the same effect becomes apparent but the profiles are no longer symmetric.

Nonlinear density profile
Now that we have verified the finite difference method we can consider different density strata. In particular, we are interested in the behaviour as the linear gradient transitions into a two-layer strata with a region of lighter fluid sitting on top of a region of heavier fluid. In order to consider the case of general stratification of the form ρ(ζ ) = f (ζ ), we use the density profile considered by Farrow and Hocking [26] in nondimensional form as Fig. 7 The nonlinear density profiles for different densities d 2 = 0.4π (left) and d 2 = 0.49π (right). This is a typical density profile for a stratified fluid as used by Farrow and Hocking [26] where the scaling has been to make the density vary by a single unit from top to bottom so that the density at the base is ρ = 1 and at the top ρ = 0 and the nondimensional density perturbation is one unit. The values of d 1 and d 2 determine the shape of the density profile, as seen in Fig. 7, and are related by the expression d 2 = arctan(d 1 ).
After differentiating ρ with respect to Ψ , we obtain Substituting (18) into (6), we find We now treat (19) as a boundary value problem with the boundary conditions (5) plus dΨ dξ → 0 as ξ → ∞ and 0 < ζ < 1, which dictates the flow becomes horizontal, and solve using the finite difference method.
For lower G, discharge into the sink is very high, while for higher G, discharge into the sink is low. Figure 8 shows streamlines for d 2 = 0.4π in which the density stratification is close to linear (see Fig. 7a), for G = 2 and G = 8 and the behaviour is very much like that shown in Fig. 2 as might be expected. Figure 9 shows the discharge profiles for the same cases and again is very similar to the linear gradient case. In the slower flow of Fig. 9b, the top layers are moving quite quickly towards the point sink while the bottom region is almost stagnant. Figure 10 clearly shows the effect of the different density stratification. In (a), d 2 = 0.4π is an almost linear stratification while in (b), d 2 = 0.49π is almost two distinct layers (see Fig. 7b). Discharge profiles are shown for the case G=1, which is a relatively strong flow, and clearly shows that the density layers in the step stratification are restricting the flow. Discharge above the interface is much stronger and this is also reflected in the streamlines in Fig. 8, where the greater separation near the base is indicative of the slower flow.
As G increases and the discharge decreases, the stratification begins to have a greater effect. Figure 11 shows the density contours about the interface with d 2 = 0.4π . In (a), G=2 and the flow is quite strong and so fluid is withdrawn over the full height, simply pulling the interface up into the outlet. However, in Fig. 11b, where G = 8, the flow is quite weak and it turns out that the upstream density gradient adjusts to having the interface much higher. Attempts to impose the condition that the interface remained at the mid-height resulted in spurious oscillations in the solutions even when the domain of the computations was greatly increased. This is an indication that there is no solution with the centre of the density interface at the middle height, thus violating the upstream condition. The reason for this is that the flow suction is not strong enough to overcome the potential energy of such a stratification.  Fig. 7. Since the density is constant along a streamline, the colours are also indicative of the density. However, the density values fit the functions shown in Fig. 7 rather than being evenly spaced (see Fig. 11) It is most likely that flows with that strength and an upstream interface at the mid-level will become single-layer flow as described later in this work. However, solutions here assume that the streamlines all enter the sink directly. The greater separation of the streamlines as ξ → ∞ is a reflection of the fact that the flows near the bottom are slowing. The case in which the bottom layer becomes separate is considered in the next section where a dividing streamline forms leaving a stagnant region near the bottom.
In cases where the stratification is more strongly layered, that is for d 2 > 0.4π say, there is a distinct layering in the flow velocity, see Fig. 10. Therefore, if the flow rate is sufficiently small, then there is no longer a solution in which the initially prescribed density structure exists upstream, and the numerical scheme adjusts to a new upstream condition in which the centre of the interface is no longer located at ζ = 0.5.
Given the adjustment of the upstream level of the interface, it is appropriate to define G * = G (1 − H mid ), where G * represents the value of G when the length scale adjusted to be the distance to the middle of the interface, and H mid is defined as the location of the density isoline with nondimensional density ρ = 0.5. H mid also coincides with the point at which the density gradient is steepest. This value will provide an estimate of the effective maximum G value. Figure 12 shows G versus G * for several different values of d 2 . Clearly as G increases, G * levels off,   Fig. 12 G vs. G * for d 2 = 0.45π, 0.4π, 0.35π , bottom to top. As G increases, G * levels off at around G * ≈ 2, and as d 2 increases the maximum G * decreases slightly, i.e. a sharper interface results in a lower minimum value of G * indicating that for any density stratification and sink height, there is a flow rate beneath which the withdrawal pressure is not large enough to pull the water beneath the interface upward into the outlet. The result of this is that for G * greater than this value, the water will almost exclusively come from this layer adjacent to the sink. In this example, G * approaches G * = 2 as G increases. This will be investigated further later in this work.

Changing sink elevation-nonlinear density profile
Thus far we have only considered situations with nonlinear stratification for which the point sink is situated on the upper surface of the region. However, it is also of interest to consider cases in which the sink is situated at some arbitrary height. The flow will still be axisymmetric, and so the only change required to the method is to modify the values of ψ on ξ = 0 so that boundary conditions are given in (7) along with Ψ ξ = 0 as ξ → ∞ to ensure the streamlines remain horizontal. Figure 13 shows the streamlines for the case with h s = 0.5, i.e. mid-height of the region, for two cases of flow rates, G = 2 and G = 8, all with d 2 = 0.4π . The relatively broad interface in this example is centred at ζ = 0.5, so that the withdrawal point is in the middle of this interface. At the higher value of G (lower flow rate), the narrowing of the flowing velocity profile can be seen. Figure 14 shows the velocity profiles for these two cases. At low flow rates, G = 8, the profiles are very narrow in the vertical direction as the interface acts like a very strongly stratified region, so the withdrawal layer is very narrow with a relatively high velocity centre. The higher flow, G = 2, looks more like the earlier results as the flow is much stronger over the whole height of the channel and the strength of the flow overcomes the stratification.
This effect is greatly exaggerated in Fig. 15, in which the density interface is the case for d 2 = 0.45π and hence is much sharper. In this example, the withdrawal layer is narrow even in the higher flow, G = 2 case, but for G = 8, the low flow case, the layer is very thin with a very sharp looking velocity profile.
If the sink height, h s , is moved away from the mid-height of the density interface, H mid , then the flow begins to separate more strongly into the two different layers: above the interface and below the interface. Figure 16b shows the case with the sink located at h s = 0.7, i.e. near the middle of the upper layer, and a sharper interface with d 2 = 0.49π . It is clear that the flow coming from below the interface is very low and for h s = 0.7, there is a clear separation in velocity between the two layers, and within each layer, the velocity profiles are quite flat, indicating almost a potential flow. In the case of the withdrawal from the middle of the interface at h s = 0.5, the flow is evenly spread over the whole region.   These results again show that the nature of the stratification has a profound impact on the level from which the fluid is drawn. The presence of an interface at any level will limit the extent of withdrawal from a part of the flow domain, and therefore needs to be kept in mind during extraction of water from aquifers, or oil from reservoirs.
When the sink is not on the top, i.e. h s = 1, H mid will eventually rise up to the level of h s as G increases, and so an appropriate calculation for the maximum G for each h s is to compute when the middle of the interface lifts off from the height ζ = 0.5, i.e. the point at which the solution is no longer valid as a flow with the upstream level of the middle of the interface at ζ = 0.5.
The correct boundary condition upstream for flows with the interface centred around h = 0.5 is that the mid-level density line and the equivalent streamline remain at this height as ξ → ∞. As G increases, this condition is violated as the flow slows. In other words, the solution is modified to one in which the interface is no longer centred at ζ = 0.5 but at some higher level. An examination of solutions with truncation of the region at larger values of ξ indicated that the streamlines do level off and the error is not due to the finite domain. Attempting to enforce the correct boundary condition resulted in the method failing to converge. The conclusion is that solutions with this upstream condition do not exist for larger values of G, i.e. slower flows. Figure 17 shows the maximum G for each h s at different d 2 . These values can be considered as a maximum for each sink height. It is interesting that the value gets lower as the interface becomes sharper. This figure is important for comparison with the single-layer flows. It is likely that at values of G larger than this, a region develops that is stagnant and not drawn towards the sink.

Axisymmetric withdrawal with flow in only one layer
In the first part of this paper, we have considered the situation of a variable density gradient and determined that if the flow rate is not large enough then solutions with flow across the full height of the channel do not exist. Therefore, we now seek solutions for steady withdrawal through a point sink in the upper layer of a two-layer fluid with an extremely thin interface in which the lower layer is stagnant. The goal is to determine the flow rate at which coning [17] occurs, in which the interface is drawn up directly into the point sink. It is to be expected that this may coincide Fig. 18 Axisymmetric flow through a horizontally confined aquifer between z = 0 and z = 1 on the free surface z = η(r ) with a sink at z = h s with the point at which the strongly layered solutions in the previous section fail as G increases. A sketch of such a flow is shown in Fig. 18, where the distance from the free interface to the top of the region is D and the flow into the sink is Q as before and H * s is the sink height above the interface. Define piezometric head functions φ 0 = p ρ 0 g +ẑ and φ 1 = p ρ 1 g +ẑ for the upper and lower layers in coordinates (r ,ẑ). Matching the pressure across the interface and nondimensionalising the φ values with respect to g D, where g = (ρ 1 − ρ 0 )/ρ 0 )g, and the lengths with respect to D gives us similar parameters to the earlier sections. Implementing Darcy's law (2) and conservation of mass for the (nondimensional) functions Φ 0 and Φ 1 in the flowing layer in nondimensional coordinates (r, z), leads to a requirement that throughout the regions above and below the interface, respectively. If the lower layer is stagnant, then Φ 1 = 0 throughout, so that our only concern is the value of Φ 0 in the upper layer, subject to conditions on the upper boundary and the interface between the two layers. Consequently, we drop the subscript and use Φ henceforth to represent Φ 0 . The nondimensional variables (r, z) are slightly different to those in the earlier section, (ξ, ζ ) because in the earlier section the whole aquifer height was scaled to be one, while here it is the distance from the top to the density interface as r → ∞. The free surface is defined as z = η(r ), and the boundary and surface conditions are where the first two conditions are applied on the free interface to ensure constant pressure in the stagnant lower fluid and that there is no flow through the surface, while the third stipulates there can be no flow through the impermeable upper surface. The function Φ must satisfy the condition that where M is representative of the sink strength. This places a mathematical point sink at (r, z) = (0, h * s ), where h * s = H * s /D, to model the withdrawal behaviour. The '*' indicates the slight change in definition of sink height to be from the level of the interface.
The quantity M is related to G in Eq. (10) by the relation where d is the total height of the confined aquifer as defined in Fig. 1, and in general will be d ≈ 2D depending on the level of the interface as r → ∞.
Care must be taken in comparing results because in the single-layer flow, all of the flux remains above the interface, while in the earlier sections, the mid-point of the interface enters the sink directly, so that only half of the flux into the point sink is from above the interface.
To find a solution to this single-layer problem, an appropriate choice for Φ that satisfies (20) and (23) is where J 0 (θr ) is a Bessel function of the first kind, and C(θ ) is a weighting function to be determined so that Φ will satisfy the surface conditions, and defines a point sink at height, z = h * s , and also satisfies (23). The integral in Eq. (26) is not amenable to a numerical solution of the nonlinear problem, and so we choose to replace it with a sum, so that we represent where θ k , k = 0, 1, 2, . . . are appropriately chosen eigenvalues.
To ensure that the interface becomes horizontal as r → ∞ in the discrete form (28), it is necessary that Φ z = 0 at r = L and so θ k are the zeros of J 0 (θ k L) = 0, for k = 1, 2, 3, . . ., and L is a large value of r at which the domain must be truncated.

Linear solution
If the flow rate is comparatively small, then (22) can be approximated with on z = 0 and also from (21) Then noting and applying the linear condition (29) in (28), gives To compute C k , multiply (32) by J 0 (θ k r )r dr and integrate from 0 to L and thus by orthogonality of the eigenfunctions (see Abramowitz and Stegun [24]), Once the C k s are known, we can use (28) Fig. 19 for the case h * s = 0.5 and with M = 0.05, 0.1 and 0.15. Experiments with the numerical calculations showed that the value of truncation L had almost no influence on the solution once L was greater than L ≈ 15, and this was also found to be the case in the nonlinear solutions below. It is possible to extend this method to solve the full nonlinear problem for the shape of the interface. Choosing the potential Φ to take the form given in (26), we must now implement the surface conditions directly on z = η(r ) rather than on z = 0. As above, an appropriate choice for Φ is given by (28).
Since z = η(r ) is unknown, we require an expression that is consistent with Φ in (28), so we let and therefore where the θ k , k = 0, 1, 2, 3, . . . are defined as above. To solve the nonlinear problem, we take these forms for Φ and η and substitute into Eqs. (21) and (22). Note that (23) is satisfied automatically by the choice of the functions for Φ.
Truncating the series to N terms gives 2N equations for the B k and C k , k = 0, 1, . . . , N − 1. Integration was performed using highly accurate Gaussian quadrature. Using the iteration routine fsolve in Octave [25], the values of the coefficients can be found.
Some examples of interface shapes are given in Fig. 19, which shows a comparison between the approximate linear solutions and the full nonlinear solutions. It is clear that for small values of flow rate M agreement is very good, but as M increases, the two solutions begin to differ as we would expect. However, the excellent comparison of the two for small values of M is confirmation of the numerical scheme.
As the value of M increases, the surface pulls up towards the sink point. At some point, coning will occur, meaning that the interface will pull directly into the sink (draw up). Figure 20 shows the height of the middle of  That two-layer solutions with a stagnant lower layer exist for values of flow rate M below M max is consistent with the existence of solutions flowing over the full depth for G 2 only below some value in the stratified fluid, i.e. above some high flow rate. However, calculations indicate that the qualitative values are not particularly close, with the value of G computed as the minimum of the single-layer flow being quite a bit larger than that obtained as the maximum of the stratified computations, suggesting further work is needed to consider the transition of the smoothly stratified case to the two-layer case.

Summary
We have considered axisymmetric flow of a stratified fluid into a point sink in a vertically confined aquifer. First we considered the case where the density stratification was linear. Two methods of solution were used, an eigenfunction expansion in Whittaker functions and then a finite difference method. A comparison showed good agreement between the two methods, verifying the numerical approach.
The finite difference method was extended to consider nonlinear density strata and different sink locations. Velocity profiles were plotted to understand the effects of variations in flow rate and density. The differences in flow behaviour for nonlinear density profiles were identified. Regions in parameter space were found in which there were no steady solutions with streamlines filling the entire solution domain. In particular, low flow rates with strongly nonlinear density gradients seemed to be limited in the scope of possible solutions. It would seem likely that there are solutions in which the streamfunction Ψ = 0 does not go along the base of the region, but instead separates from the wall with a stagnant region below. However, this calculation would be difficult using this finite difference method. Another method would be necessary to consider the formation of this stagnant region in situations with a general density stratification.
This suggestion that there is such a stagnant region led us to consider the case of flow in only a single layer above a sharp interface, and we were able to show not only that such a layer exists, but also identify the maximum flow rates beneath which they can be found. The next step is to identify such solutions with a variable density strata.

Funding Information 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/.