Selective withdrawal of a three-layer viscous fluid

Two dimensional, selective withdrawal of a three-layer fluid in a rectangular tank is considered. Each fluid is viscous, weakly compressible and miscible. The lowermost, heaviest fluid flows with constant speed out through line sinks located on the base and lower side walls of the tank. The uppermost and lightest fluid recharges the tank via line sources located along the top and on the upper side walls at a rate that matches the outward flux. If the line sinks are turned on suddenly, the interfaces between the fluids are drawn downward at different rates; this disparity arises due to the fact that the wave configurations that propagate along the two interfaces are not the same. Density and vorticity contours are presented and used to determine the time at which the middle transition layer of fluid begins to be withdrawn from the tank; at this point the rate of extraction of the most dense fluid is its greatest. Inviscid linearised interfaces are calculated which serve as checks on the early stage evolution of the full nonlinear model.


Introduction
The extraction of fluid from a rectangular reservoir is a problem that has been well studied in recent decades due to its many engineering applications. For example, Imberger and Hamblin [1] explore the different mechanisms that may contribute to setting the observed state of a body of water; in their example the concern was with predicting the location of the interface between two different density fluids when a river flows into a lake. Dufresne et al. [2] described the flow patterns that may arise in shallow rectangular reservoirs which are often used in the design of storm water management systems. The applications are not limited to water flow examples, as Cohen et al. [3] illustrated how the extraction of a fluid from a finite reservoir may be optimised as a mechanism to coat micro-particles. This process is important for local drug delivery and is used in a plethora of medical settings.
The overwhelming majority of literature discussing selective withdrawal focuses on configurations in which two fluids of different densities are separated by a single, sharp interface. It is widely recognised that this is a simplification of many real-life situations; rarely are physical problems described by two fluids separated by a sharp interface. In many circumstances the two fluids blend more gradually from one to the other and this is the situation that motivates the work described here. As a simple model of this type of problem we consider selective withdrawal in a three-layer fluid where the middle layer serves to play the role of a buffer between the dense lower fluid and the relatively light upper fluid.
A key question pertaining to withdrawal problems was posed by Lister [4]. Suppose that a body consists of fluid whose density varies with height (or depth) and that a sink is introduced that extracts fluid from the system. Lister [4] asked whether the withdrawn material has a density equal to that of the fluid immediately adjacent to the sink, or whether the extraction process is in itself sufficient to deform the density surfaces so that the removed fluid is drawn from a mixture of different densities? Lister proved that the presence of surface tension is necessary to guarantee the stability of the interface when only a single fluid layer is being extracted. Additionally, computations enabled a critical flow rate to be predicted and, at least for super-critical flow, it was possible to determine the proportion of each fluid extracted from the system. Forbes and Hocking [5] considered a rectangular tank containing two fluid layers of different densities with the more dense fluid extracted from the base of the tank by means of a line sink(s). They showed, at least for ideal fluids drawn through a single line sink, that at early times the position of the interface between the two fluids can be predicted using a linearised theory but as the interface nears the base of the tank a nonlinear model is required. Moreover Farrow and Hocking [6] showed that the drawdown of an interface between two fluids can occur over a range of flow rates and proved that the finite length of the tank is an important factor that can either delay or promote the fluid extraction. They also found that if the line sink sits on the base of the tank and is turned on impulsively, activated waves begin to propagate along the interface. These waves reflect from the side walls of the tank and increasing the depth of the bottom layer of fluid allows the transient wave activity to decay before the drawdown commences.
Cosgrove and Forbes [7] built on the studies of Farrow and Hocking [6], Forbes and Hocking [5], Hocking [8] and others, but extended the analysis to allow for weakly compressible, viscous and miscible fluids that flow rotationally. The motivation for this work was the suggestion that such non-ideal fluids might align more closely with real world situations. The aim of Cosgrove and Forbes was to withdraw as much of the more dense, lower fluid as possible before the interface between the two fluids reached the line sink(s). Both symmetrically and asymmetrically placed sinks were considered and it was proven that two symmetrically placed line sinks on the base of the tank maximised the withdrawal of the lower fluid. Moreover, it was shown that adding drains to this configuration tended to reduce the formation of waves at the interface of the two fluids, thereby increasing the available time to withdrawal the lower fluid before the interface reached the drains.
The selective withdrawal of a viscous fluid into a line sink on the side wall of a rectangular tank was described by Imberger, Thompson and Fandry [9]. They explain that the reflection of shear waves off the rear wall can significantly modify the flow and that the length of the channel and the Froude number are both of critical importance. Furthermore at large Prandtl numbers the layers tend to become unsteady and collapse occurs quickly. It is suggested that in a real world situation the base of the tank or reservoir is unlikely to be flat and probably would be sloped or of irregular shape; this would naturally further complicate the wave structure although they provided no data to quantify this.
Çalışkan and Elçi [10] investigated selective withdrawal from reservoirs in an attempt to assess water quality. They considered a three-dimensional hydrodynamic model, as this best represents a real reservoir. Withdrawal from different depths within the body of water was modelled and it was concluded that placing the line sink at the bottom of the reservoir, as opposed to suspending it somewhere within the fluid, was the preferred choice to optimise the quality of the extracted water.
Information regarding the selective withdrawal of a three-layered fluid from a rectangular tank is rather limited, and what is known tends to concern various solitary wave solutions. These nonlinear solutions can only be realised on large bodies of fluid for which edge effects are either small or are simply ignored. Such solitary waves cannot arise in a tank of finite width as used in this paper. Much of the motivation for examining three layer systems arises in the context of modelling of ocean or atmospheric fluids. For example, Molin et al. [11] investigated the 'sloshing' of a rectangular tank containing three fluids which has a focus on the mixing of the boundary layers when the tank is subject to different forced motions such as swaying and rolling. Davey [12] considered a three-layer fluid model and examined its susceptibility to baroclinic instability which typically is seen in stably stratified rotating fluids [13].
The three-layer intrusion flow model considered by Chen and Forbes [14] might be regarded as occurring in an unbounded rectangular tank containing a middle fluid layer sandwiched between two thin interfaces in the vertical direction. All three horizontal fluid layers are in motion with propagating waves fixing the shape of the two interfaces. They used a perturbation expansion strategy to develop a linearised solution associated with small wave amplitudes and implemented spectral methods to compute several nonlinear solutions for steady periodic waves. In this paper we employ similar techniques to generate both linear and nonlinear solutions. We consider the withdrawal of fluid which is initially at rest and set in motion by the impulsive opening of three line sinks. One of these is located on the base of the rectangular tank while the other two are positioned at the bottom of each side wall. Although different drain configurations are considered, the problem closely follows solutions presented in [7], [6] and [5] with the key modification that we have three fluid layers of distinct densities which are separated by two sharp interfaces. The fluids are assumed to be weakly compressible, viscous and flow rotationally. Spectral decomposition is used to identify nonlinear solutions pertaining to various combinations of tank geometries and flow rates. A linearised solution is obtained which serves as a benchmark to confirm the validity of the non-linear approximation at early times.
The remainder of this paper is organised as follows. The mathematical model which describes the problem is formulated in Sect. 2, in which equations for the inviscid model are developed along with a linearised solution of the problem in Sect. 2.1 and the spectral method that is used to perform nonlinear calculations is described in Sect. 2.2. The results of these computations are set out in Sect. 3 in the form of density and vorticity contour plots. Finally, in Sect. 4 we close with a few comments and some concluding remarks. The chosen geometry, shown in Fig. 1, was largely motivated by the findings of Cosgrove and Forbes [7]. In their work they concluded that two line sinks placed symmetrically on the base of the tank was the optimal positioning for the withdrawal of the lower, denser fluid layer to be maximised. This prompted us to place line sinks on the lower side walls of our tank together with a centrally located line sink on the base of the tank; our aim was to further enhance the withdrawal of the lower fluid layer by 'flattening out' the lower interface. One unfortunate consequence of this choice is that it leads to an ill-posed Neumann problem. Pivato [15] states the circumstances under wihch Laplace's equation can be solved when Neumann conditions are prescribed on all the boundaries of a region. In simple terms, this requires the line integral of the flux constraint to vanish on each of the boundaries and this can be ensured by the introduction of suitable line sources at the top of each side wall, as described above.
At time t = 0 all three drains are opened simultaneously and fluid leaves at a constant rate of 2V 0 from the sinks located at the bottom of the side walls and at a rate V 0 from the sink on the base of the tank. To ensure the container remains full at all times, light fluid is admitted to the tank via the uppermost line source at a constant rate of AV 0 /W and through the side-wall sources at a rate 2V 0 . This choice of flow rates ensures that the inward and outward fluxes balance (each have a value 4 AV 0 ). Assuming that all three fluids are weakly compressible, fluid on either side of the interfaces undergoes a smooth, but rapid, transition from one density regime to the next. We therefore take the density in the form written ρ(x, y, t) = ρ 2 +ρ(x, y, t) where ρ 2 is taken as the reference density andρ(x, y, t) is a small perturbation. Furthermore, we assume that all three fluids have a viscosity ν and diffuse at a rate σ .
We are now in a position to express the problem in a convenient dimensionless form. We write all lengths relative to the depth of the lowermost layer H , which throws up the natural scales √ H/g for the time and √ g H for fluid speeds; here g denotes the acceleration due to gravity. Densities are scaled on that of the middle fluid ρ 2 . The problem is therefore reduced to one that depends on the eight parameters (2.1) A sketch of the dimensionless tank is depicted in Fig. 1. The parameter λ measures the depth of the middle fluid relative to the other two layers so that the overall height of the tank is 2 + λ and initially the two interfaces are located at y = 1 and y = 1 + λ. The width of the tank is 2β and the sink at the centre of the base is of extent 2α while the sources and sinks on the sidewalls are all of width α/2. A measure of the flow rate is provided by the Froude number F while the two coefficients D 1 (< 1) and D 2 (< 1) describe the density ratio between fluids that share an interface. Finally, the parameters μ and γ are the dimensionless viscosity and diffusion coefficients respectively.

The inviscid model
Following the arguments given by Forbes and Hocking [5] and Farrow and Hocking [6], we proceed under the assumption that each of the fluids flow irrotationally. The velocity vector can be written as the gradient of a suitable velocity potential q = u k i˜+ v k j = ∇ φ k for k = 1, 2, 3; here u and v denote the horizontal and vertical velocity components of the flow in the direction of the i˜and j unit vectors, respectively.
At times t > 0 we shall suppose the upper and lower interfaces to be positioned at y = η 1 (x, t) and y = η 2 (x, t) respectively. They move and deform owing to the inward flow from the line sources at the top of the tank and the outward flow from the line sinks beneath. It follows that Laplace's equation is satisfied in each of the three regions, and therefore The presence of line sources on the top and line sinks on the bottom of the side walls, require that the following boundary conditions are satisfied at x = ±β: On the base of the tank y = 0 we require that while at the top y = 2 + λ the recharge rate is given by As the fluid within each region is unable to cross an interface into an adjacent zone, it is necessary to impose suitable kinematic conditions at each of the interfaces. Moreover, the fluid pressure must be continuous across each of the interfaces, see [5], which leads to the two dynamic conditions It is also necessary to allow for time-dependent velocity potentials written in the form where j is the (yet to be determined) time-dependent flow which is superimposed on the time-independent φ b background flow. This base field is determined by examining the flow in the tank should it contain a single constant density, non-viscous, incompressible fluid that moves irrotationally. This sets the steady state behaviour in the region −β < x < β, 0 < y < 2 + λ and is given by in which the coefficients are given by (2.10) It is easily verified that the time-dependent potentials j (x, y, t) satisfy Eqs. (2.2) with φ j is replaced by j for j = 1, 2, 3. Moreover the boundary condition on the top of the tank (2.5) becomes the boundary condition on the bottom of the tank, (2.4) is while the boundary conditions on the side walls of the tank (2.3) are simply The given problem reduces to one of determining the two time-dependent interface functions η 1 (x, t) and η 2 (x, t) together with the three time-dependent velocity potentials j (x, y, t) ( j = 1, 2, 3). The highly nonlinear nature of the problem means that full results rely on the implementation of a suitable numerical method but a linearised solution is possible. This serves twin purposes; on the one hand it gives a guide as to the nature of the solution at small flow rates and, on the other, it provides a description of the early stage evolution of fully nonlinear structure. We consider this linearised solution next.
Forbes and Hocking [5] obtained a linearised inviscid solution for a tank that contained two fluids separated by a single interface. We remark that although we shall consider nonlinear simulations in later sections of the paper, these inviscid linearised considerations are helpful markers that assist in verifying the early evolution of the full nonlinear problem. The presence of a third fluid and the second interface introduces some complications which is encapsulated within a coupled system of equations. We begin by expanding the velocity potentials in terms of the (assumed small) Froude number F so that (2.14) where φ j L represents the first order velocity potentials for j = 1, 2, 3. In a similar vein the two interface functions can be expressed to first order as The leading order potentials satisfy with the associated boundary conditions derived from (2.11), (2.12) and (2.13). This gives (2.18) The reduced velocity potentials assume the general forms By taking the appropriate partial derivatives of (2.14) and making use of the two interface functions (2.15) and (2.16), the linearised dynamic conditions on the upper and lower interface are represented by The substitution of the appropriate time derivatives of expressions (2.19) into the linearised interface requirements (2.20) and (2.21) leads to (2.23) We then use (2.15) and (2.16) in conjunction with expansions of the horizontal and vertical velocity components u k and v k taken either side of the two interfaces. The substitution of the results into the kinematic conditions (2.6) generates two expressions for each of the linearised interface functions η 1L and η 2L . If Eqs. (2.22) are differentiated with respect to time and compared to the appropriate expressions derived from the kinematic equations, four relations are obtained. Two of these relate to the upper interface so that while the two for the lower interface give in which n1 (x) = a n sinh cos nπ x β , n2 (x) = a n sinh cos nπ x β , a n = (2.26) We can now deduce the key equations for the functions K n (t) and L n (t). We remark that expressions (2.23) involve the time derivatives of Q n (t), R n (t), S n (t) and T n (t) and these quantities can be found by differentiating equations (2.24) and (2.25) with respect to time and solving for the four unknowns. This yields and the substitution of (2.27) into (2.23) leads to the desired coupled system of second order, constant coefficient ordinary differential equations. These may be written as where The system (2.28) can be rewritten as a system of four first order equations K n (t) = X 1 (t), L n (t) = X 2 (t), subject to the initial conditions where 1n and 2n are given by (2.26). If we march the solution of this system forward in time we obtain the forms of K n (t) and L n (t) which can be used in conjunction with Eqs. (2.22) to find the terms in the expressions for the two interface positions correct to order F as given by (2.15) and (2.16).

The viscous model
In our formulation of the linearised model solutions were obtained by solving Laplace's equation in each of the three fluid regions. When extending to the nonlinear model, we adopt the strategy outlined by Cosgrove and Forbes [7] and derive solutions by considering the entire tank to be a single fluid region with an initial prescribed density profile. With the two-dimensional flow q = ui˜+ v j , where u and v denote the horizontal and vertical velocity components, Boussinesq theory allows us to express the conservation of mass in terms of an incompressible part plus an effective transport condition ∂ρ ∂t + q · ∇ρ = γ ∇ 2ρ , (2.30) in whichρ ≡ρ(x, y, t) denotes a small perturbation to the density ρ and γ is the appropriate dimensionless diffusion parameter. The streamfunction ≡ (x, y, t) can now be introduced. It too can be considered as the sum of two components; one being the steady state background flow b ≡ b (x, y) to which we add an unsteady viscous correction c ≡ c (x, y, t). We establish the background part by ignoring the effects of viscosity and considering the steady flow of a single fluid in the region |x| <= β, in which the time-independent Fourier coefficients A n and B n are given by (2.10). Following Farrow and Hocking [6] and Cosgrove and Forbes [7], slip conditions are imposed along the side walls and on the base of the tank so that the boundary layers on these surfaces can be ignored. As the boundary conditions (2.3), (2.4) and (2.5) are satisfied by the steady component b , the unsteady correction term has to account for the viscosity of the fluid and we write where the time dependent Fourier coefficients A mn (t) are to be determined. Therefore, the full representation of the streamfunction takes the form (2.31) As described in [7], the horizontal and vertical velocity components of the viscous momentum equations can be expressed in terms of vorticity as a single scalar equation, given by in which the vorticity is represented as and μ is the dimensionless viscosity parameter. As the streamfunction identically satisfies the incompressible component of the mass equation ( where the time dependent Fourier coefficients, B mn (t), are currently unknown. The arbitrary function of integration f (y, t) was written as a Fourier series which is even-valued about y = 0 so that this particular form was chosen in order to minimise the amount of Gibbs' phenomenon present in the solutions. Therefore, just as described in Sect. 2, the full representation of the density function takes the form The important evolution equations for the time-dependent Fourier coefficients A mn (t) and B mn (t) are derived using spectral decomposition. The resultant system of equations are similar to those presented in Cosgrove and Forbes [7]; and for ease of reference for the reader, we include our results in Appendix A.
The fluid is initially assumed to be at rest so that all the coefficients A mn (0) = 0. The starting values for the coefficients B mn follow from the form of the initial density profile given bȳ Once we had deduced the various Fourier coefficients it was then a simple task to reconstruct the forms of the streamfunction, vorticity and density using (2.31), (2.36) and (2.38) respectively.

Results and discussion
In all our calculations we set the diffusion and viscosity parameters μ = γ = 10 −4 for consistency with the values used in [7]. While the qualitative behaviour of the results is unlikely to be sensitive to the values of γ and μ, they would need to be adjusted should there be a particular physical problem of interest. For example, should there be temperature or salt diffusion, it would be appropriate to set γ ≈ 0.1μ or γ μ, respectively. Extensive testing suggested that the number of Fourier coefficients M = N = 31 provided a good compromise between the need to generate accurate solutions while ensuring that computational run times were not excessive. To check the adequacy of the resolution we set M = N = 51 in many of our calculation combinations and in all cases there was no discernible change to the results which suggests that our computations were sufficiently resolved; we generate all the figures with this number of Fourier coefficients. We supposed that the three fluid densities are given by ρ 1 = 0.900, ρ 2 = 0.950 and ρ 3 = 1.000 so that the density ratios across the upper and lower interfaces are D 1 = ρ 1 /ρ 2 = 0.9474 and D 2 = ρ 2 /ρ 3 = 0.9500 respectively. Informed by the methods used by Cosgrove and Forbes [7] we optimised the choices of the various parameters by varying each in turn in order to achieve maximal total outflow V . For fixed α the total outflow flux can be estimated as V = 4α Ft w , where t w is the time taken for the lower interface to reach one or more of the line sinks at the base of the tank. Furthermore, we do not have a precise definition for the interfaces because the three distinct density layers are modelled by a smeared density profile across the depth of the tank. For computational purposes we estimated the location of each interface to be Table 1 Total outflow V and corresponding withdrawal time t shown for various combinations of α and F. Other parameters are fixed with β = 8 and λ = 1 the position where the calculated density equals the average of fluid densities either side of that interface; in other words the locations where ρ = 0.925 and 0.975 were used to track the upper and lower interfaces respectively. Figure 2 depicts the relationship between the total outflow V and the Froude number F for various values of α. It appears that the outflow is optimal when F ≈ 0.03; a result that seems to hold over the entire range of α values tested. The justification for this choice of F is given by the observation that decreasing F much below 0.03 led to significantly longer computation times with an almost negligible increase in the outflow. At the other extreme, with F > 0.08, there was virtually no outflow. Table 1 lists the withdrawal times and total outflow for various combinations of α and F. It is remarked that the value of α = 1 appears to be close to the optimal drain half-width in the sense that this choice seems to achieve high outflow in a relatively short time. This tentative result can be justified by comparing the various outputs corresponding to F = 0.03 (the previously determined optimal Froude number). Although Table 1 shows that the total outflow V is greater for α = 0.95 than α = 1, this marginally increased outflow of only 1% is achieved at the expense of increasing the withdrawal time, t w , by more than 6%. Conversely, if α is increased to 1.05, the withdrawal time, t w , is reduced by 11% in comparison to α = 1, but this time saving is accompanied by a reduction in the total outflow V of greater than 6%. Thus, based on these observations, α = 1 provides a good balance between computation time and outflow.
Density contour plots were generated and some representative solutions are shown in Fig. 3. Here the contours show the variation in densities in the vicinity of the interfaces. The white space in the following density plots denotes regions in which the fluid density equals the undisturbed values corresponding to the upper, middle and lower fluid regions. In all cases it was found that the lower interface was drawn through the lower side drains before reaching the central drain on the base of the tank. When α = 1 the extraction times t w were determined for both F = 0.03, 0.07 (see Table 1) by determining when the the lower interface began to be drawn through the uppermost parts of the sidewall drains, as shown in Fig. 3. This is a good visual confirmation that increasing the Froude number much beyond the optimal F = 0.03 leads to a more rapid extraction of fluid and with a concomitant reduction in the total outflow. We remark that the small zone of high density contours ρ = 0.975 that appear towards the top of Fig. 3b, are spurious artifacts caused by the numerical approximation, probably as a result of the Gibbs' phenomenon. As such they can be safely ignored.
Having decided that the optimal Froude number is in the vicinity of F = 0.03 and with α = λ = 1,we next looked at the effect of varying the tank half-width β. Some representative results are shown in Fig. 4 which shows that the withdrawal time t w is well approximated via the linear relationship t w = (β −0.51)/0.09. Figure 5 shows the extraction times for the two values β = 2 and β = 6. We point out that these results hold for a fixed dimensionless height of the tank, so care must be taken if these findings are extended to examine other geometrical combinations.
Remembering that as the draining started impulsively, it might be expected that both the upper and lower interfaces would be drawn toward the line sinks located on the bottom of the side walls and on the base of the tank. Figure 6 shows density contours plotted using (2.38) with the inviscid linearised interfaces (2.15) and (2.16) superimposed. Both interfaces are initially drawn toward the sinks located at the bottom of the walls, mainly owing to the inward linearised interfaces lose accuracy but this in unsurprising as nonlinear effects increasingly become the dominant mechanism. Shown in Fig. 7a are contours of vorticity, plotted using (2.36) with the optimal parameter values α = 1 and F = 0.03 and with β = 8 and λ = 1. Figure 7 shows vorticity contours immediately after the line sources and sinks are activated, where small but relatively intense regions of vorticity can be observed ; the production of vorticity is likely due to the bending of the isopycnal contours, as described by the Dubreil-Jacotin-Long equation [16]. In the central region there are small relatively weak vorticity structures forming on the lower interface, indicating that fluid is beginning to be drawn out through the line sink on the base of the tank. As time progresses through Figs. 7b-d the regions of vorticity on the sides of the tank expand in size but decrease in magnitude, in contrast, the vorticity centrally located on the lower interface expand and deepen. This is reassuring confirmation of the behaviour suggested by the density profiles in Fig. 6. Just as for the density plots, the linearised interfaces are in good accord with the vorticity behaviour at early times but diverge as time progresses and nonlinearity becomes more prominent.
One salient question raised in our study relates to the importance of modelling selective withdrawal with a three-layer rather than a two-layer profile. In order to assess the possible importance of the additional layer, we performed a few experiments in which the density of the upper most layer was set equal to that of the middle layer, 6 Linearised interfaces (dashed lines) superimposed on the computed density contour plots at times t = 3, 6, 12 and 24. The tank parameters are chosen to be α = 1, β = 8, λ = 1 and F = 0.03 namely ρ 1 = ρ 2 = 0.95, effectively simplifying the problem to that of two fluids separated by a single interface. We simulated the two-layer system using the parameter values α = 1, β = 8, λ = 1 and F = 0.03 and the results are illustrated in Fig. 8. With these parameter choices, Fig. 8a shows the interface of the two-layer fluid achieves a maximum withdrawal time of t w = 54, with a corresponding total outflow of V = 6.48. At this time, Fig. 8b shows the lower interface of the three-layer fluid is still well above the line sinks on the lower side walls so this solution can be continued; in fact it may be taken on until t w = 82 and produces a total outflow of V = 9.84 as noted in Table 1. Therefore, for this particular case, the three-layer model produces a fluid output in excess of one-third greater than that of the equivalent two-layer counterpart. To ensure these results were not in some way peculiar and atypical, we simulated several other cases taken from Table 1 and repeated the calculations when the densities of the upper two layers were set equal. In all cases the outflow V was markedly greater in the case of the three-layer fluid system. We conclude that, in general, a three-layer fluid system exhibits significantly different behaviour to that of its two-layer counterpart. It is also worth noting that the linearised interfaces for our chosen tank geometry exhibit vastly different behaviour from that of [5,7], most notably our linearised interfaces appear to remain almost 'flat'. Changes in the value of λ affect both the depth of the middle fluid layer together with the height of the tank. Figure 9 shows that as the value of λ falls so the top interface begins to assume the shape of the lower interface, as seen in Fig. 9d. This indicates that as λ → 0 the three-layer problem approaches that of a two-layer case as we almost have two fluids separated by a single interface. However, this hypothesis could not be confirmed by taking λ to very small values, because then proper resolution of the thin middle layer requires increasing numerical resolution.

Conclusion
In this paper we have considered the selective withdrawal of a three-layer viscous fluid from a two-dimensional rectangular tank, with the aim of estimating the maximum volume of the most dense fluid that can be drawn from the tank before the undesirable middle transition layer begins to be removed as well. A suite of calculations gives the empirical result that the optimal flow rate occurs when the Froude number is approximately F = 0.03 and the drain half-width α = 1. At small Froude number a linearised solution can be developed by treating the pair interfaces as a coupled system; it is comforting that this gives a good description of the flow, at least for early times. After an impulsive start to the draining process, both linearised interfaces were drawn toward line sinks at the bottom of the side wall, with the lower interface more noticeably drawn toward the sink on the base of the tank initially. Travelling waves then propagate along the interfaces towards the tank walls, effectively reversing the drawdown effects. Over time the linearised interface continues to be drawn towards the base of the tank, but it descends at a much slower rate than seen in the corresponding nonlinear solutions; this would be anticipated on physical grounds. We point out that this is in stark contrast to the behaviour observed in two-layer problems [5,7], when there seems to be no significant difference in the speeds at which the interface moves. We also remark that our calculations comparing two and three layer configurations suggest that inclusion of the third layer can have a surprisingly large effect on withdrawal times and outflow rates.
Nonlinear solutions were obtained using spectral methods. Contour and vorticity plots were employed as visualisation tools for studying the maximum withdrawal of the most dense, lower layer of fluid; the nonlinear predictions were pleasingly in accord with the linearised results at early times. In comparison to analogous two-layer problems, the interfaces of the three-layer system remained quasi-horizontal while drawn down occurred, with the lower interface being drawn out of the lower side wall drains before reaching the line sink on the base of the tank. This behaviour suggests that the addition of a transition layer into a selective withdrawal model might lead to improved quality of extracted fluids which may be significant for various industrial applications such as storm water management [2]. Moreover it was observed that as the Froude number increases so the interfaces are more rapidly drawn toward the line sinks on the lower side walls.
Much of what has been described here provides a foundation for further studies into the withdrawal of liquid from three-layer configurations. While, as expected, some of the properties of the flow are reminiscent of two-layer results, there are important extra features that do not appear to have been noted in the earlier studies. It would be of interest to investigate the robustness of some of these novel structures and ask whether they are inherent three-layer effects, or present only in particular flow regimes. Furthermore, here we have deliberately restricted ourselves to one particular geometrical configuration; can better withdrawal properties result from moving the sinks and sources to more advantageous locations? These, and other related issues, are deferred and await further work.  in which the constants C kl are defined by (2.37). This system is advanced in time using a standard fourth order Runge-Kutta method starting from prescribed initial conditions.