1. Introduction
The Hele-Shaw flow or the gap-averaged Stokes flow is an important subclass of fluid problems, where the flow of fluids occurs between two closely placed plates. In such a case, one ignores the out-of-plane velocity component and averages the in-plane velocity over the thickness of the gap to reduce the problem to two dimensions. The Hele-Shaw flow attracts considerable attention because of its applications in oil recovery (Hornof & Baig Reference Hornof and Baig1995; She et al. Reference She, Aoki, Wang, Li, Nasir, Mahardika, Patmonoaji, Matsushita and Suekane2022), micro-fluidics (Hashimoto et al. Reference Hashimoto, Garstecki, Stone and Whitesides2008; Chakraborty et al. Reference Chakraborty, Ricouvier, Yazhgur, Tabeling and Leshansky2019) and porous media flow (Saffman & Taylor Reference Saffman and Taylor1958; Taylor & Saffman Reference Taylor and Saffman1959), etc. In the oil recovery process, for example, one might consider that the oil is getting extracted through a sink while being surrounded by air or a different less viscous fluid that tries to penetrate the oil and therefore has an impeding effect on the recovery process. The flow problem is quite challenging – one has to track one or multiple moving interfaces that typically exhibit Saffman–Taylor-like instability (Saffman & Taylor Reference Saffman and Taylor1958). As the system evolves, the interfaces develop viscous fingers giving rise to beautiful and complex patterns (Paterson Reference Paterson1981; Chen Reference Chen1989; Li et al. Reference Li, Lowengrub, Fontana and Palffy-Muhoray2009; Zhao et al. Reference Zhao, Ying, Lowengrub and Li2017, Reference Zhao, Li, Ying, Belmonte, Lowengrub and Li2018).
The classic Hele-Shaw flow has been studied extensively through experimental (Paterson Reference Paterson1981; Chen Reference Chen1989), numerical (Li, Lowengrub & Leo Reference Li, Lowengrub and Leo2007; Zhao et al. Reference Zhao, Ying, Lowengrub and Li2017; Morrow, Moroney & McCue Reference Morrow, Moroney and McCue2019; Morrow, De Cock & McCue Reference Morrow, De Cock and McCue2023) and analytical (Escher & Simonett Reference Escher and Simonett1996, Reference Escher and Simonett1997; Prokert Reference Prokert1998; Tanveer & Xie Reference Tanveer and Xie2003; Xie & Tanveer Reference Xie and Tanveer2003) means. Extensions of the classical problem have also been formulated and investigated, where the nature of the fluid (Kondic, Palffy-Muhoray & Shelley Reference Kondic, Palffy-Muhoray and Shelley1996; Fast et al. Reference Fast, Kondic, Shelley and Palffy-Muhoray2001), the geometry of the system (Dias & Miranda Reference Dias and Miranda2013) and driving forces (Miranda & Widom Reference Miranda and Widom2000; Zhao et al. Reference Zhao, Niroobakhsh, Lowengrub and Li2021, Reference Zhao, Anjos, Lowengrub, Ying and Li2023; Anjos et al. Reference Anjos, Zhao, Lowengrub and Li2022) have been varied. The literature is extensive, and we do not intend to give a comprehensive review here. In the current work, our interest is in observing such flows in radial cells but in a multi-layer set-up with a sink as the driving force. We discuss a few key references below.
Various experimental, numerical and analytical studies have been conducted for the multi-layer Hele-Shaw problem. For example, the annular flow was considered experimentally (Cardoso & Woods Reference Cardoso and Woods1995) and in a rotating radial cell (Carrillo, Soriano & Ortın Reference Carrillo, Soriano and Ortın1999, Reference Carrillo, Soriano and Ortın2000). It was found that any perturbation to the outer interface tends to stabilize as the interfaces approach one another and the annulus region gets thinner (Cardoso & Woods Reference Cardoso and Woods1995). In a separate work, an empirical relation between capillary number and another dimensionless quantity (related to the ratio of centrifugal to capillary forces) was found for a wide range of values (Carrillo et al. Reference Carrillo, Soriano and Ortın1999). A linear stability analysis, carried out in conjunction with experiments, revealed a good match between theoretical and experimental observations for the number of fingers produced on the interface (Carrillo et al. Reference Carrillo, Soriano and Ortın2000). In an early study, the onset of Rayleigh–Bénard convection in presence of magnetic fields was checked (Aniss, Brancher & Souhar Reference Aniss, Brancher and Souhar1993). Logvinov (Reference Logvinov2019) investigated the displacement of a more viscous fluid with a less viscous one in presence of a source. Through linear stability analysis, the author identified a mode that grew the fastest. Also, the predictions matched quite well with the experimental results in the low capillary number regime.
In the analytical front, the use of complex variable theory has proven quite fruitful as the real and complex parts of analytic functions are harmonic. Taking a cue from this, powerful techniques have been devised (Richardson Reference Richardson1996; Crowdy Reference Crowdy1999, Reference Crowdy2002; Crowdy & Kang Reference Crowdy and Kang2001; Cummings & King Reference Cummings and King2004). The effect of surface tension is ignored in these studies as it is not easy to find solutions in the presence of the capillary forces. More recently, attention was given to the annular problems using a pressure differential (Dallaston & McCue Reference Dallaston and McCue2012). Again, surface tension was ignored to bring in the force of complex analysis. In contrast, a number of other works consider the multi-layer Hele-Shaw problem with surface tension (Beeson-Jones & Woods Reference Beeson-Jones and Woods2015; Gin & Daripa Reference Gin and Daripa2015, Reference Gin and Daripa2021; Anjos & Li Reference Anjos and Li2020). Instead of using the complex variable approach, all of them took up a small-perturbation analysis approach to investigate the stability of the interfaces. For example, authors tried to find the optimal value of the viscosity of the intermediate fluid in order to inject fluid at the fastest rate possible while not disrupting the stable flow (Beeson-Jones & Woods Reference Beeson-Jones and Woods2015). In a separate work, an upper bound on the growth of perturbations was found and verified with simulations (Gin & Daripa Reference Gin and Daripa2015). The scope of this analysis was expanded further with analysis carried out for a three layer Hele-Shaw problem with the middle layer having variable viscosity (Gin & Daripa Reference Gin and Daripa2021). The goal was to find injection schemes that would maintain the stability of the interfaces. Prior to that, the question of short-time existence and uniqueness of the Hele-Shaw problem for various initial conditions of the interface and in the presence of surface tension was settled (Escher & Simonett Reference Escher and Simonett1996). Anjos & Li (Reference Anjos and Li2020) used a second-order mode coupling theory to demonstrate that as the thickness of the annulus domain decreases, the interaction between the interfaces gets strong, with wide fingers forming on the interfaces with bifurcated tips. However, they observed that if the thickness of the annulus is reduced too much, then the finger-splitting morphologies are replaced by polygonal-like structures with narrow fingers.
As mentioned above, most of the analytical studies either disregard the effects of the surface regularization mechanism, i.e. capillary effects, or prove the existence and uniqueness of the solutions for a short time, or rely on a perturbation approach that linearizes the problem and can therefore again be relied upon for a short duration of time. In some recent analytical works (see Green, Lustri & McCue Reference Green, Lustri and McCue2017), attention has been given to surface tension; however, the geometric set-up has been special in those works, requiring certain symmetries. This is where the role of the numerical methods becomes important. For example, in a very recent work, Morrow et al. (Reference Morrow, De Cock and McCue2023) performed simulations in annular domains with surface tension using level set methods to understand the progression of viscous fingering. The motion was driven by either rotation or pressure difference. A sharp interface approach involving the boundary integral technique yields much accurate results for a long time duration of the problem, especially when the space–time convergence of the problem is of high order. For example, Zhao et al. (Reference Zhao, Anjos, Lowengrub and Li2020) have investigated the pattern formation problem for a three-layer problem driven by a source at the origin using the boundary integral approach. Nonlinear simulations are shown to match with experiments and weakly nonlinear analysis, although the simulations go well beyond the weakly nonlinear regime and provide good insight to fully nonlinear dynamics. The question of numerical stiffness in the time stepping, due to interfacial conditions, is dealt with using a small-scale decomposition technique (Hou, Lowengrub & Shelley Reference Hou, Lowengrub and Shelley1994). There are several advantages of using the boundary integral method. In a moving boundary problem such as ours, the method allows us to recast the original problem formulated in terms of partial differential equations to boundary integral equations defined only on the interfaces, and then to track the latter. This leads to a reduction of the problem dimension by one. Other benefits include the exact treatment of the interface conditions and the existence of highly accurate numerical techniques. A major downside of the boundary integral method is its inability to deal with topological changes in the configuration.
An early numerical work in the presence of an eccentric sink is by Kelly & Hinch (Reference Kelly and Hinch1997), who consider the problem in the presence of small surface tension. Using a boundary integral method, they show that a zero surface tension cusp formation scenario can be avoided even if the surface tension effect is small. The scope of this work is expanded further in Tian & Nie (Reference Tian and Nie1998) and Ceniceros, Hou & Si (Reference Ceniceros, Hou and Si1999). Both use boundary integral formulations like Kelly & Hinch (Reference Kelly and Hinch1997); however, the numerical methods differ. One of the main contributions of Tian & Nie (Reference Tian and Nie1998) is to analyse the nature of the singularity in a sink driven flow. They predict that the interface reaches the sink before all fluid is sucked out, which is also supported by the experimental evidence found in Paterson (Reference Paterson1981). In Ceniceros et al. (Reference Ceniceros, Hou and Si1999), the authors compute cases of (a) small and (b) large viscosity ratios between the outer and inner fluids. It is observed that when viscosity ratio is small, the interface develops a finger that evolves into a wedge having a neck region. In the case where the ratio is large, the formation of the neck gets suppressed, and the finger that develops on the interface is thinner.
The current work is motivated by the cusp-like interface morphology that develops during sink-driven Hele-Shaw flow (Tian & Nie Reference Tian and Nie1998). We are interested in investigating the effects of the sink on the dynamics of the Hele-Shaw problem with one more interface, thus the sink location could be in the interior of the inner interface or the annulus region, as shown in figure 1. It should be noted that in the latter case, one does not have any results from linear analysis. This makes the problem difficult to solve via analytical means, hence in this work, we adopt a numerical approach based on boundary integral formulations (Zhao et al. Reference Zhao, Anjos, Lowengrub and Li2020). To the best of our knowledge, this is the first boundary integral theory based work applied to the multi-connected Hele Shaw problem driven by a sink. Through numerical simulations, we observe that the interaction between the two interfaces introduces rich dynamics beyond the classic cusp-like patterns for a single interface (Tian & Nie Reference Tian and Nie1998). Our study reveals the importance of initial distance between the two interfaces in the nature of pattern formation – if the two interfaces are initially placed ‘close’, then they tend to come close to each other before either one of them reaches the sink; however, if they are ‘well separated’ at the beginning, then one interface reaches the sink before the other catches up with it. This leads to two typical events: (1) a cusp-like pattern forming mechanism if one of the two interfaces reaches the sink faster than the other; (2) interface-merging patterns if they come close to each other. In particular, we observe that the inner interface can be wrapped by the other to have both scenarios. We find that multiple parameters contribute to the dynamics, including the width of the annular region, the location of the sink and the mobilities of the fluids. An important practical application of the current study could be in the oil extraction process where multiple layers of oil get recovered through the sink with air or water trapped in the oil. The success of the process would depend on how the viscous fingers form. For example, if the fingers generated by the air bubble reach the sink before the oil is extracted, then the recovery efficiency might be reduced. Two-interface problems driven by a Darcy-type equation in a geometric set-up like ours can also found in other areas, such as tumour dynamics (Lu et al. Reference Lu, Hao, Liu, Lowengrub and Li2022), where, moving from inner to outer region, we find necrotic core, tumour and healthy tissues, respectively. Thus the current problem connects to other areas of application and is of fundamental importance.
Schematic diagram of a three-layer Hele-Shaw flow in the presence of a sink. The innermost layer is 
$\varOmega _1$, which is surrounded by a domain 
$\varOmega _2$ of annulus-like shape. The outermost layer is 
$\varOmega _3$. The moving interface between 
$\varOmega _1$ and 
$\varOmega _2$ is 
$\varGamma _1(t)$, and that between 
$\varOmega _2$ and 
$\varOmega _3$ is 
$\varGamma _2(t)$. The sink is represented by the black dot. We always place the sink at the origin. Depending on the location of 
$\varOmega _1$, the sink can be either (a) in the fluid region 
$\varOmega _1$ or (b) in the annulus 
$\varOmega _2$.
The paper is organized as follows. In § 2, we describe the governing equations. In § 3, we discuss our numerical methods. In § 4, we discuss the main results. In § 5, we summarize our findings.
2. Governing equations
We consider a radial Hele-Shaw cell with three fluid layers trapped between two plates separated by a small distance 
$b$, which remains unchanged. The innermost fluid region 
$\varOmega _1$ is a bounded, simply connected domain in 
$\mathbb {R}^2$. The region 
$\varOmega _1$ is surrounded by a second fluid that occupies an annulus-like region 
$\varOmega _2$, and 
$\varOmega _2$ in turn is surrounded by a third fluid domain 
$\varOmega _3$, which extends to infinity. The closed interface that separates 
$\varOmega _1$ and 
$\varOmega _2$ is denoted by 
$\varGamma _1(t)$, and the one that separates 
$\varOmega _2$ and 
$\varOmega _3$ is denoted by 
$\varGamma _2(t)$, as shown in figure 1.
In each of these regions, the fluid is considered to be incompressible and irrotational. Therefore the gap-averaged velocity follows Darcy's law
\begin{equation} {\boldsymbol{u}}_{i} =-\frac{b^2}{12 \mu_i}\,\boldsymbol{\nabla} P_i =-\bar{M}_i\, \boldsymbol{\nabla} P_i, \quad {\boldsymbol{x}} \in \varOmega_i, \end{equation}
where 
${\boldsymbol {u}}_{i}$ is the velocity, 
$P_i$ is the gap-averaged pressure, 
$\mu _i$ is the viscosity of the fluid, and 
$\bar {M}_i$ is the mobility in the domain 
$\varOmega _i$, 
$i=1,2, 3$. The incompressibility condition requires
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}}_{i} = 0, \quad {\boldsymbol{x}} \in \varOmega_i. \end{equation}
The irrotational nature of the velocity fields, i.e. 
$\boldsymbol {\nabla } \times {\boldsymbol {u}}_i=0$ in the fluid domains, implies that the problem can be recast in terms of a velocity potential in each of the domains that satisfies the Laplace equation there.
In the present problem, the flow is driven by the removal of the fluids through the sink following the equation
\begin{equation} -Q = \int_{\varSigma_0} {\boldsymbol{u}}_k\boldsymbol{\cdot} {\boldsymbol{n}} \,{\rm d}s, \end{equation}
where 
$Q$ is the net flux out of the system. Here, 
$k$ is either 1 or 2, corresponding to the location of the sink in 
$\varOmega _1$ or 
$\varOmega _2$. Also, 
$\varSigma _0$ is a small interface around the point of extraction, mimicking the existence of a tube that is used to extract the fluids, 
${\boldsymbol {n}}$ is the unit outward normal on 
$\varSigma _0$, and 
$s$ is the arc length of the interface. Note that the point of extraction is always at the origin of the system. Depending on our goal, we suitably adjust the geometry to place the extraction point in either 
$\varOmega _1$ or 
$\varOmega _2$.
The pressure is discontinuous across the two interfaces:
\begin{equation} P_2-P_3 = \sigma_{23} \kappa_{23}\ \text{on}\ \varGamma_2(t),\quad \text{and}\quad P_1-P_2 = \sigma_{12} \kappa_{12}\ \text{on}\ \varGamma_1(t), \end{equation}
where 
$\sigma _{12}$ and 
$\sigma _{23}$ are the surface tensions, and 
$\kappa _{12}$ and 
$\kappa _{23}$ are the curvatures of the interfaces 
$\varGamma _1(t)$ and 
$\varGamma _2(t)$, respectively. The kinematic conditions or the continuity of the normal components of the fluid velocities on the interfaces read
$$\begin{gather} {\boldsymbol{u}}_2\boldsymbol{\cdot} {\boldsymbol{n}}= {\boldsymbol{u}}_3\boldsymbol{\cdot} {\boldsymbol{n}}\quad \text{on } \varGamma_2(t), \end{gather}$$
$$\begin{gather}{\boldsymbol{u}}_1\boldsymbol{\cdot} {\boldsymbol{n}}= {\boldsymbol{u}}_2\boldsymbol{\cdot} {\boldsymbol{n}}\quad \text{on } \varGamma_1(t). \end{gather}$$
We use the length scale 
$L_0=R_1(0)$ (initial size of the inner interface) and the time scale 
$T_0=2 {\rm \pi}\,R_1^2 (0)/Q$ to non-dimensionalize the system. We obtain the non-dimensional equations
$$\begin{gather} \boldsymbol{u}_i=-M_i\,\boldsymbol{\nabla} P_i, \quad \text{for}\ \boldsymbol{x}\in \varOmega_{i}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_i=0,\quad \text{for}\ \boldsymbol{x}\in \varOmega_{i}, \end{gather}$$
$$\begin{gather}P_1-P_2=\frac{1}{Ca}\,\kappa_{12}, \quad \text{for}\ \boldsymbol{x}\in\varGamma_1(t), \end{gather}$$
$$\begin{gather}P_2-P_3=\frac{\alpha}{Ca}\,\kappa_{23}, \quad \text{for}\ \boldsymbol{x}\in\varGamma_2(t), \end{gather}$$
$$\begin{gather}\boldsymbol{u}_1\boldsymbol{\cdot}\boldsymbol{n}=\boldsymbol{u}_2\boldsymbol{\cdot}\boldsymbol{n}, \quad \text{for}\ \boldsymbol{x}\in\varGamma_1(t), \end{gather}$$
$$\begin{gather}\boldsymbol{u}_2\boldsymbol{\cdot}\boldsymbol{n}=\boldsymbol{u}_3\boldsymbol{\cdot}\boldsymbol{n}, \quad \text{for}\ \boldsymbol{x}\in\varGamma_2(t), \end{gather}$$
$$\begin{gather}{\int_{\varSigma_0}\boldsymbol{u}_k\boldsymbol{\cdot} {\boldsymbol{n}}\,{\rm d}s}=-2{\rm \pi},\quad k=1 \text{ or } 2, \end{gather}$$
where the capillary number 
$Ca$ indicates the relative importance of the viscous to surface tension forces:
\begin{equation} Ca=\frac{ Q\,R_1(0)}{2 {\rm \pi}\sigma_{12} M_0}. \end{equation}
Here, 
$M_0$ is a characteristic mobility, and 
$M_i =\bar {M}_i /M_0$ is the dimensionless mobility of the 
$i$th fluid. The parameter 
$\alpha = \sigma _{23} / \sigma _{12}$ is the ratio of the surface tensions. We define a new function 
$\phi _i = -M_i P_i$ such that 
$\Delta \phi _i =0$ in each fluid domain. We will formulate the numerical method using this function 
$\phi _i$. Equation (2.13) is the scaled version of (2.3). Note that in the non-dimensionalization, we scale out the strength 
$Q$ of the sink. Finally, an experimental work with real fluids used the ratio of the surface tensions of the interfaces 
$\alpha = 0.485$ (Cardoso & Woods Reference Cardoso and Woods1995). For simplicity, we assume that 
$\alpha =1$ in this paper, though we can use different surface tension parameters 
$\sigma _{23}$ and 
$\sigma _{12}$ in our simulations.
3. Boundary integral formulation and time-stepping algorithm
Since the function 
$P_i$ (or 
$\phi _i = -M_i P_i$) is harmonic in 
$\varOmega _i$, using the potential theory, we rewrite the boundary value problem in terms of integrals:
\begin{equation} \phi_i(\boldsymbol{x})=\frac{1}{2{\rm \pi}}\int_{\varGamma_1}\gamma_1\,\frac{\partial \ln |\boldsymbol{x}-\boldsymbol{x}'|}{\partial \boldsymbol{n}(\boldsymbol{x}')}\,{\rm d} s(\boldsymbol{x}')+\frac{1}{2{\rm \pi}}\int_{\varGamma_2}\gamma_2\,\frac{\partial \ln |\boldsymbol{x}-\boldsymbol{x}'|}{\partial \boldsymbol{n}(\boldsymbol{x}')}\,{\rm d} s(\boldsymbol{x}')-\ln|\boldsymbol{x}|,\end{equation}
where the first two terms correspond to the double-layer representation of a harmonic function 
$\phi _i$ in the the fluid domain 
$\varOmega _i$, using two unknown dipole densities 
$\gamma _1$ and 
$\gamma _2$. The density functions 
$\gamma _1$ and 
$\gamma _2$ are defined on the boundaries 
$\varGamma _1(t)$ and 
$\varGamma _2(t)$, respectively. The effect of a sink has been incorporated in the solution by the term 
$-\ln |{\boldsymbol {x}}|$ (Greenbaum, Greengard & McFadden Reference Greenbaum, Greengard and McFadden1993; Zhao et al. Reference Zhao, Anjos, Lowengrub and Li2020). In two dimensions, the Green's function 
$G({\boldsymbol {x}}, {\boldsymbol {0}} ) = -\ln |{\boldsymbol {x}}|$ is harmonic in 
$\mathbb {R}^2\setminus \{{\boldsymbol {0}}\}$ and satisfies the equation 
$-\Delta G = 2{\rm \pi} \, \delta ({\boldsymbol {x}})$, where 
$\delta ({\boldsymbol {x}})$ is the Dirac delta function at the origin (also the location of the sink).
Using the pressure jump conditions across the interface, we obtain a system of integral equations for the unknown density functions 
$\gamma _1$ and 
$\gamma _2$:
$$\begin{gather}
\frac{1}{2}\left(\frac{M_2}{M_1}+1\right)\gamma_1(\boldsymbol{x}_1)+\frac{1}{2{\rm \pi}}\left(\frac{M_2}{M_1}-1\right)
\left( \int_{\varGamma_1
}\gamma_1(\boldsymbol{x}')\,\frac{\partial \ln
|\boldsymbol{x}_1-\boldsymbol{x}'|}{\partial\boldsymbol{n}(\boldsymbol{x}')} \,{\rm d}
s(\boldsymbol{x}')\right.\nonumber\\ +\left.
\int_{ \varGamma_2
}\gamma_2(\boldsymbol{x}')\,\frac{\partial \ln
|\boldsymbol{x}_1-\boldsymbol{x}'|}{\partial
\boldsymbol{n}(\boldsymbol{x}')}\,{\rm d}
s(\boldsymbol{x}')-2{\rm \pi} \ln
|\boldsymbol{x}_1|\right) =-\frac{1}{Ca}\,\kappa_{12},
\end{gather}$$
$$\begin{gather}\frac{1}{2}\left(\frac{M_2}{M_3}+1\right)\gamma_2(\boldsymbol{x}_2)+\frac{1}{2{\rm \pi}}\left(1-\frac{M_2}{M_3}\right)\left(
\int_{\varGamma_1
}\gamma_1\,(\boldsymbol{x}')\frac{\partial \ln
|\boldsymbol{x}_2-\boldsymbol{x}'|}{\partial
\boldsymbol{n}(\boldsymbol{x}')}\,{\rm d}
s(\boldsymbol{x}')\right.\nonumber\\ + \left.
\int_{\varGamma_2
}\gamma_2(\boldsymbol{x}')\,\frac{\partial \ln
|\boldsymbol{x}_2-\boldsymbol{x}'|}{\partial
\boldsymbol{n}(\boldsymbol{x}')}\,{\rm d}
s(\boldsymbol{x}')- 2{\rm \pi}\ln
|\boldsymbol{x}_2|\right) =-\frac{1}{ Ca}\,\kappa_{23}.
\end{gather}$$ These two equations are Fredholm integral equations of the second kind, well-conditioned from a computational point of view. Both the integral operators in (3.2) and (3.3) are compact, and the kernels have a removable singularity. Once the integral equations are solved and dipoles 
$\gamma _1$ and 
$\gamma _2$ are obtained, one can use the Dirichlet–Neumann map to compute the normal velocities of the interfaces as
\begin{align} V_{\varGamma_1}&=\frac{1}{2{\rm \pi}}\int_{\varGamma_1}\gamma_{1,s'}\,\frac{(\boldsymbol{x}-\boldsymbol{x}')^\perp\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})}{|\boldsymbol{x}-\boldsymbol{x}'|^2}\,{\rm d}s'(\boldsymbol{x}') \nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{\varGamma_2}\gamma_{2,s'}\,\frac{(\boldsymbol{x}-\boldsymbol{x}')^\perp\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})}{|\boldsymbol{x}-\boldsymbol{x}'|^2}\,{\rm d}s'(\boldsymbol{x}')- \frac{\boldsymbol{x}\boldsymbol{\cdot}\boldsymbol{n}}{|\boldsymbol{x}|^2}, \end{align}
\begin{align} V_{\varGamma_2}&=\frac{1}{2{\rm \pi}}\int_{\varGamma_1}\gamma_{1,s'}\,\frac{(\boldsymbol{x}-\boldsymbol{x}')^\perp\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})}{|\boldsymbol{x}-\boldsymbol{x}'|^2}\,{\rm d}s'(\boldsymbol{x}') \nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{\varGamma_2}\gamma_{2,s'}\,\frac{(\boldsymbol{x}-\boldsymbol{x}')^\perp\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})}{|\boldsymbol{x}-\boldsymbol{x}'|^2}\,{\rm d}s'(\boldsymbol{x}')-\frac{\boldsymbol{x}\boldsymbol{\cdot}\boldsymbol{n}}{|\boldsymbol{x}|^2}, \end{align}
where the subscript 
$s$ denotes the partial derivative with respect to arc length, and 
${\boldsymbol {x}}^{\perp }=(x_2,-x_1)$. The interfaces evolve through these velocities.
The integral equations are solved following a Nyström method whereby the integral equations are discretized at marker points 
${\boldsymbol {x}}_i$ using spectrally accurate quadrature rules. Since the kernels of integral equations are periodic and smooth, the trapezoidal rule with modified kernels has spectral convergence. One can also use an alternating point quadrature rule to achieve the same effect (Sidi & Israeli Reference Sidi and Israeli1988). The resulting linear system is solved via the generalized minimal residual (GMRES) method (Saad & Schultz Reference Saad and Schultz1986). The integral operators in the Dirichlet–Neumann maps can similarly be computed with the same accuracy, making the overall numerical computation spectrally accurate in space. A core component of GMRES requires computing the matrix–vector product. Since in our case the matrix is dense but structured, one can use the fast multipole method (Greengard & Rokhlin Reference Greengard and Rokhlin1987) or fast tree-code (Lindsay & Krasny Reference Lindsay and Krasny2001; Feng et al. Reference Feng, Barua, Li and Li2014) to expedite the computation. This reduces the cost of matrix–vector products from 
$\textit {O}(N^2)$ to 
$\textit {O}(N\log N)$ or even 
$\textit {O}(N )$, where 
$N$ is the size of the matrix (the total number of marker points).
One fundamental challenge in the surface-tension-driven Hele-Shaw flow is how to update the interface efficiently and accurately. A straightforward analysis of the equations of motion shows that one has to maintain the condition 
$\Delta t \sim \Delta x^3$ if the time-stepping method is explicit. Here, 
$\Delta t$ and 
$\Delta x$ are the sizes of the time step and space resolution, respectively. Satisfying this stability constraint requires very small time steps. The computational cost gets really high, especially for complicated interfaces where a large number of marker points are needed to maintain good spatial resolution.
The small-scale decomposition (SSD) technique (Hou et al. Reference Hou, Lowengrub and Shelley1994) alleviates the problem and reduces the spatio-temporal constraint to 
$\Delta t \sim \Delta x$. Following this technique, we first rewrite the equation of motion in terms of the length 
$L(t)$ of the interface and the tangent angle 
$\theta _j$ that the marker point 
${\boldsymbol {x}}_j$ makes with the positive direction of the 
$x$-axis. The equation for arc length 
$L$ is non-stiff and can be updated using the second-order Adams–Bashforth method. However, the 
$\theta$ equation is stiff. Following the idea of SSD, we recast the equation in terms of a stiff part and a non-stiff part. The stiff part is linear and can be integrated exactly in the Fourier space. For completeness, we briefly describe the method here.
The algorithm requires the marker points to be equally spaced in arc length at all times. This is achieved through a direct discretization at 
$t=0$ and by adding a special tangential velocity 
$T_{\varGamma _i}$, 
$i=1,2$, of the form
\begin{equation} T_{\varGamma_i}(\alpha,t) = T_{\varGamma_i}(0,t)+ \int_{0}^{\alpha} s_{\alpha^{\prime}}\kappa V_{\varGamma_i} \,{\rm d}\alpha^{\prime} - \frac{\alpha}{2{\rm \pi}}\int_{0}^{\alpha} s_{\alpha^{\prime}}\kappa V_{\varGamma_i} \,{\rm d}\alpha^{\prime} \end{equation}
to the equations of motion at later times to maintain the equal space property, where 
$\alpha$ parametrizes the interface, and 
$s_\alpha =\sqrt {x_{\alpha }^2+y_{\alpha }^2}.$ Also, 
$V_{\varGamma _i}(\boldsymbol {x}(\alpha,t))$ and 
$T_{\varGamma _i}(\boldsymbol {x}(\alpha,t))$ denote the normal velocity and tangential velocity of the interface 
$\varGamma _i(t)$.
In the 
$({\boldsymbol {s}},{\boldsymbol {n}})$ (tangent–normal) frame, the equations of motion then become
\begin{equation} \frac{{\rm d}}{{\rm d}t} {\boldsymbol{x}}=V_{\varGamma_i}(\boldsymbol{x})\, \boldsymbol{n}+T_{\varGamma_i}(\boldsymbol{x})\,\boldsymbol{s}, \quad \textrm{for}\ \boldsymbol{x}\in\varGamma_i(t), \ i=1,2, \end{equation}
where 
$\boldsymbol {n}$ and 
$\boldsymbol {s}$ represent the unit normal and tangential vectors on each interface. Next, using the equal arc length frame, we repose the equations of motion in terms of 
$L$ and 
$\theta$ coordinates as
$$\begin{gather} \theta_t = V_s + \kappa T, \end{gather}$$
$$\begin{gather}s_{\alpha,t} = \left(T_s - \kappa V \right)s_{\alpha}, \end{gather}$$
where we use the relation 
${\rm d}s = ({L}/{2{\rm \pi} })\,{\rm d}\alpha$. In both of these equations, we suppress 
$\varGamma _i$ to keep the notation simple. Also, 
$t$ in the subscript of different variables denotes derivative with respect to time. The first equation is stiff, while the second equation is not and can be updated using an explicit scheme, e.g. the second-order Adams–Bashforth method. Following (Hou et al. Reference Hou, Lowengrub and Shelley1994), we recast (3.8) in the form
\begin{equation} \theta_t = \frac{\sigma}{s_{\alpha}^3}\,\mathcal{H}[\theta_{\alpha \alpha \alpha} ] + N(\alpha, t), \end{equation}
where the operator 
$\mathcal {H}[\,\cdot\, ]$ denotes the Hilbert transform, and 
$N(\alpha, t) = V_s + \kappa T - (\sigma /s_{\alpha }^3)\,\mathcal {H}[\theta _{\alpha \alpha \alpha } ]$ is non-stiff and has a removable singularity. In the Fourier space, it can be diagonalized as
\begin{equation} \hat{\theta}_t(t) =-\frac{\sigma k^3}{s_{\alpha}^3}\,\hat{\theta}(k,t) + \hat{N}(k,t). \end{equation}
We implement a linear-propagator-based Adams–Bashforth scheme of second-order accuracy to numerically integrate this equation, then perform an inverse Fourier transform to find 
$\theta$. We also use smoothing filters and cut-off filters to control the onsets of non-physical high-frequency spurious modes (Jou, Leo & Lowengrub Reference Jou, Leo and Lowengrub1997).
4. Results and discussions
In the following subsections, we investigate various mechanisms of instability in the three-layered configuration. The two interfaces move due to removal of fluid through a sink, located in either the fluid domain 
$\varOmega _1$ or 
$\varOmega _2$. Unless stated otherwise, we use 
$N=8192$ and 
$\Delta t=1\times 10^{-5}$, where 
$N$ is the number of marker points on each interface, and 
${\Delta } t$ is the time step. The iterative GMRES solver tolerance is set to 
$\epsilon = 10^{-12}$, so is the filter tolerance 
$\epsilon _{{tol}}$.
All our computations are carried out using a computer with 3.7 GHz AMD Ryzen ThreadRipper 3970X CPUs. We start our simulations using smooth interfaces with relatively simple geometries. At the beginning, it takes only a few GMRES iterations to obtain the solution. However, all our simulations approach finite time singularities, and near the breakdown, the count of iterations increases dramatically due to (i) the tiny distance separating the two interfaces, (ii) the high curvature development (sharp corners), or (iii) the thin neck formation of the interface.
The capillary number 
$Ca$ is a dimensionless quantity representing the relative effect of viscous drag forces versus surface tension forces acting across an interface. Due to a large length scale 
$R_1(0)$ and the extraction flux time scale 
$Q$ used to non-dimensionalize the equations, our definition of the capillary number is 
$Ca={12 \mu _2 Q\,R_1(0)}/{2 {\rm \pi}\sigma _{12} b^2}$. Following this definition, for example, we find that silicone oil with viscosity 
$\mu _2=11.4$ Pa s, surface tension 
$0.02\ {\rm N}\ {\rm m}^{-1}$, gap width 
$b=0.75$ mm, initial size of the inner interface 
$R_1(0)=3$ cm, and 
$Q=0.1\ {\rm cm}^2\ {\rm s}^{-1}$ will result in 
$Ca\approx 580$ (Nase, Derks & Lindner Reference Nase, Derks and Lindner2011). One could use other less viscous fluid than the silicone oil, with smaller 
$b$ or large 
$R_1(0)$, to get this capillary number as well. In this paper, we set 
$Ca=500$ throughout.
4.1. Numerical convergence
In this subsection, we summarize the spatio-temporal convergence studies of our numerical algorithm. We introduce fluid mobility, which is widely used in porous media flow and is quite useful for further discussion. In the porous media literature, the mobility is defined as 
$M={k}/{\mu }$, where 
$\mu$ is the viscosity of the fluid, and 
$k$ is the permeability of the surrounding media. Comparing porous media and Hele-Shaw flow, we observe that the constant 
$b^2/12$ takes the role of parameter 
$k$ in the case of the latter. Hence 
$M$ varies inversely with the fluid viscosity. To study interface instabilities, we set the 
$M_{i}$ in such a way that the innermost fluid 
$\varOmega _1$ is the most viscous, followed by the annulus fluid region 
$\varOmega _2$, and the outer fluid domain 
$\varOmega _3$. In our simulations, we choose the mobilities of the fluids as 
$M_1=0.01$, 
$M_2=1$ and 
$M_3=100$ in regions 
$\varOmega _1$, 
$\varOmega _2$ and 
$\varOmega _3$, respectively. We set the initial outer interface 
$\varGamma _2(0)$ with Cartesian coordinates 
$x=({\sqrt {2}}/{6})(4\cos (\alpha )+\cos (2\alpha ))$ and 
$y=({\sqrt {2}}/{6})(4\sin (\alpha )+\sin (2\alpha ))$, where 
$\alpha \in [0,2{\rm \pi} ]$ is a parametrization. The initial inner interface 
$\varGamma _1$(0) is just a circle of radius 0.65 centred at the origin. Because of the set-up of our problem, the fluid domain 
$\varOmega _1$ gets drained from the system.
First, we demonstrate spatial accuracy. We define numerical error
\begin{equation} \text{Err}(t)=\left|A(t)-A(0)\right|, \end{equation}
where 
$A(t)$ is the area of fluid domain 
$\varOmega _2$ and should be equal to its initial value 
$A(0)$ in theory, because the location of the sink is in 
$\varOmega _1$. We plot 
$-\log _{10}\text {Err}(t)$ with respect to time 
$t$ for various values of 
$N$, the number of marker points on the interfaces. We choose 
$\Delta t=5\times 10^{-6}$. In figure 2(a), we observe that the curves are on top of each other, indicating that the solutions of the integral equation are almost identical as long as the interface is well resolved. For example, 
$N=1024$ is enough to run the simulation to 
$t=0.15$, and a further increase in the number of points does not contribute to the accuracy, suggesting the spectral accuracy of our method (Kress Reference Kress2014; Trefethen & Weideman Reference Trefethen and Weideman2014). The inset shows the region near time 
$t=0.15$ where the simulation stops. We note that the smallest distance between the two interfaces at 
$t=0.15$ is approximately 
$7.8\times 10^{-3}$, about twice the spatial resolution 
$\Delta x=4\times 10^{-3}$. Numerically, the two close interfaces result in a very ill-conditioned linear system with large condition numbers, and the GMRES iterations do not converge.
(a) Spectral accuracy of the algorithm, plotting 
$-\log _{10}\text {Err}(t)$ versus time 
$t$ for different values 
$N=1024,2048, 4096$ to see that the curves are on top of each other. The inset shows the region near time 
$t=0.15$ in detail. This is where the simulation stops because inner and outer interfaces come very close to each other. (b) Second-order accuracy of the time-stepping algorithm, plotting 
$-\log _{10}\text {Err}(t)$ versus time 
$t$ for 
$\Delta t=4\times 10^{-5}$, 
$2\times 10^{-5}$, 
$1\times 10^{-5}$, 
$5\times 10^{-6}$.
The second-order accuracy of the time-stepping scheme can also be demonstrated by using different time steps to perform the same simulation. In figure 2(b), we choose 
$N=4096$ and run four sets of simulations, with 
$\Delta t=4\times 10^{-5}$, 
$2\times 10^{-5}$, 
$1\times 10^{-5}$ and 
$5\times 10^{-6}$, i.e. each subsequent time step is half of the previous value. This suggests that when we plot 
$-\log _{10}\text {Err}(t)$ against the time 
$t$, the curves should be apart by 
$\log _{10} 4 = 0.602$ for a second-order time-stepping method, which indeed is consistent with the implemented second-order Adams–Bashforth scheme. The inset displays the final configuration of the interfaces when we stop the simulations.
4.2. Numerical validation
We consider two circles initially centred at the origin, 
$x_1^2+y_1^2=1$ (inner interface) and 
$x_2^2+y_2^2=4$ (outer interface). The mobilities of the fluids are 
$M_1=0.01$, 
$M_2=1$ and 
$M_3=100$. We choose the capillary number 
$Ca=500$, and 
$Q=-1$. For this perfect annulus problem with the sink placed right at the centre, there exist analytical solutions. Namely, the normal velocity of each interface is given by 
$V_{\varGamma _1}={1}/({\sqrt {1-2t}})$ and 
$V_{\varGamma _2}={1}/({\sqrt {4-2t}})$. In figure 3(a), we compare the numerical normal velocities from our scheme with the theoretical results, and find that they are in excellent agreement with a discrepancy of approximately 
$10^{-12}$.
(a) The normal velocity of the interfaces with respect to time. The interfaces are two circles with radii 
$1$ and 
$2$ centred at the origin. (b) The evolution of the relative perturbations 
$a_{n}(t)/R_{1}(t)$ and 
$b_{n}(t)/R_{2}(t)$. The dashed curves show the result given by the numerical approach, and the solid lines are predicted by the linear theory. The initial conditions for the inner and outer interfaces are 
$r_{1}(\alpha,0)=1.5+0.05\cos (4\alpha )$ and 
$r_{2}(\alpha,0)=2+0.1\cos (4\alpha )$, respectively. In addition, we set 
$Ca=500$, 
$M_1=0.01$, 
$M_2=1$ and 
$M_3=100$. In the case of the fully nonlinear numerical amplitudes, we utilized 
$N=8192$ points along each interface, and time step 
$\Delta t=1\times 10^{-5}$.
Next, we provide a comparison between our numerical calculations and the predictions of linear stability analysis (Beeson-Jones & Woods Reference Beeson-Jones and Woods2015; Zhao et al. Reference Zhao, Anjos, Lowengrub and Li2020; Gin & Daripa Reference Gin and Daripa2021). We consider the two interfaces as perturbed circles centred at the origin, 
$r_{1}(\alpha,t)=R_1(t)+a_n(t)\cos (n\alpha )$ (inner interface) and 
$r_{2}(\alpha,t)=R_2(t)+b_n(t)\cos (n\alpha )$ (outer interface). Here, 
$R_1(t)$ represents the size of the inner interface, and 
$a_n(t)$ denotes the cosine perturbation; 
$R_2(t)$ represents the size of the outer interface, and 
$b_n(t)$ denotes the cosine perturbation. From the linear stability analysis, the motion of the perturbations satisfies
$$\begin{gather} \dot{a}_{n} = f_{1}\left[\frac{n{\,-\,}f_{1}^{-1}}{R_{1}^2} {\,-\,} \left(\frac{M_1}{M_1{\,-\,}M_{2}}\right) \frac{n(n^2{\,-\,}1)}{Ca\,R_{1}^{3}}\right]{a}_{n} {\,+\,} f_{2}\left[\frac{n}{R_{2}^2} {\,-\,} \left(\frac{M_3}{M_{2}{\,-\,}M_3}\right) \frac{n(n^2{\,-\,}1)}{Ca\,R_{2}^{3}}\right]{b}_{n}, \end{gather}$$
$$\begin{gather}\dot{b}_{n} = f_{3} \left[\frac{n}{R_{1}^2} {\,-\,} \left(\frac{M_1}{M_1 {\,-\,} M_{2}}\right) \frac{n(n^2{\,-\,}1)}{ Ca\,R_{1}^{3}}\right]{a}_{n} {\,+\,}f_{4} \left[\frac{n{\,-\,}f_{4}^{-1}}{R_{2}^2} - \left(\frac{M_3}{M_{2}{\,-\,}M_3}\right) \frac{n(n^2{\,-\,}1)}{ Ca\,R_{2}^{3}}\right]{b}_{n}, \end{gather}$$where
\begin{equation} \left. \begin{aligned} f_{1} & = \frac{A_{12} (1-A_{23}R^{2n})}{1+A_{12}A_{23}R^{2n}} ,\quad f_{2} = \frac{A_{23} (1+A_{12})R^{(n-1)}}{1+A_{12}A_{23}R^{2n}},\\ f_{3} & = \frac{A_{12} (1-A_{23})R^{(n+1)}}{1+A_{12}A_{23}R^{2n}},\quad f_{4} = \frac{A_{23} (1+A_{12}R^{2n})}{1+A_{12}A_{23}R^{2n}}. \end{aligned} \right\} \end{equation}
Here, 
$A_{12}=({M_1-M_2})/({M_1+M_2})$ and 
$A_{23}=({M_2-M_3})/({M_2+M_3})$ are the viscosity contrasts of fluids 1 and 2, and 2 and 3, written in terms of the fluid mobilities, and 
$R=R(t)={R_1}/{R_2}$ (Zhao et al. Reference Zhao, Anjos, Lowengrub and Li2020).
For the interface configurations, we choose 
$n=4$, 
$R_1(0)=1.5$, 
$a_n(0)=0.05$, 
$R_2(0)=2$ and 
$b_n(0)=0.1$. Other parameters are 
$Ca=500$, 
$M_1=0.01$, 
$M_2=1$ and 
$M_3=100$. In figure 3(b), we plot the evolution of the relative perturbations 
${a_{n}(t)}/{R_{1}(t)}$ and 
${b_{n}(t)}/{R_{2}(t)}$ as functions of time. With the given parameters, both perturbations increase, indicating that both interfaces are unstable. The dashed curves show the results given by the numerical method, and the solid lines are predicted by the linear stability analysis in (4.2) and (4.3). The plot shows excellent agreement between the numerical and linear analysis at early times, when the perturbations are small and satisfy the assumption of linear analysis.
4.3. Motivation behind our numerical simulations
Before we discuss our results, we briefly review the important findings of a single-layer Hele-Shaw flow with suction (Tian & Nie Reference Tian and Nie1998). One starts with an initial shape of the viscous fluid domain as
\begin{equation} f(z) = \tilde{a}_1(t)\,z+\tilde{a}_2(t)\,z^2, \end{equation}
where 
$z\in \mathbb {C}$ with 
$|z|<1$, 
$\tilde {a}_1(0) = 2\sqrt {2}/3$ and 
$\tilde {a}_2(0) = \sqrt {2}/6$, and investigates the evolution under various strengths of surface tension. Then one can show that in the absence of surface tension, the interface forms a single cusp well before any part of it reaches the sink. In the cases with non-zero surface tension, the interface forms a finger that moves towards the sink. A large surface tension leads to a ‘fat’ finger, and the movement towards the sink is slow. These findings reaffirm the regularizing nature of surface tension in sink-driven Hele-Shaw flows.
Our numerical investigation starts with the interface outlined in (4.5); however, we do not restrict ourselves just to this interface. We scale up the investigation by adding the second interface and targeting its impact on the dynamics by focusing on the initial geometry of inner and outer interfaces, the location of the sink, and the effects of mobility (viscosity). We next use various configurations other than (4.5), and summarize their common characteristics in § 4.8 in terms of the evolution of surface energy. In this paper, we are interested in the interfacial instabilities. Thus we choose the outermost fluid to have the highest mobility, which makes the outer interface more unstable than the inner one. In the scheme, we can set the mobility parameters to any value. We observe that in certain cases, when we change the mobility parameters, the interfacial patterns do not change appreciably. In the next few subsections, we report some of the typical findings that we have observed.
4.4. Pattern formation with a sink and geometrically similar outer and inner interfaces
As a first variation on the classic simulations (Tian & Nie Reference Tian and Nie1998), we wish to investigate how the proximity of the outer interface to the inner one affects the pattern formation. We take the initial shape of the outer interface to be a magnified version of the inner one given by (4.5).
In figure 4(a), the initial outer interface is 1.2 times larger than the inner one, while in figure 5(a), it is 1.05 times. The mobilities of the fluids are 
$M_1 = 1$, 
$M_2 = 100$ and 
$M_3 = 10\,000$ in regions 
$\varOmega _1$, 
$\varOmega _2$ and 
$\varOmega _3$, respectively. We place the sink at the origin. In figures 4(c) and 5(c), we show our numerical results at the time when the simulation stops. The red and blue curves indicate the outer and inner interfaces, respectively.
Effects of geometry. The sink is at the origin. The mobilities of the fluids are 
$M_1=1$, 
$M_2=100$ and 
$M_3=10\,000$. In (a), we display the initial configuration. At 
$t=0$, the inner and outer interfaces are both given by (4.5). The outer interface is 1.2 times the inner interface. The intermediate and final configurations are shown in (b,c). In (d), we show a close-up of the inner interface finger at the times 
$t=0.1118$, 
$0.1128$, 
$0.1138$, 
$0.1148$ and 
$0.1150$.
Effects of geometry. The sink is at the origin. The mobilities of the fluids are 
$M_1=1$, 
$M_2=100$ and 
$M_3=10\,000$. In (a), we display the initial configuration. At 
$t=0$, the inner and outer interfaces are both given by (4.5). The outer interface is 1.05 times the inner interface. The intermediate and final configurations are shown in (b) and (c).
In figures 4(a–c), we observe the gradual formation of a finger on the inner interface that eventually approaches the sink. We have to stop the simulation at 
$T=0.1150$ due to non-convergence of the linear solver beyond this time. The part of the outer interface located near the negative 
$x$-axis and close to this finger also moves towards the sink. This result is similar to that observed in figure 2 or figure 7 of Tian & Nie (Reference Tian and Nie1998), indicating that the coupling effects of the two interfaces is weak. In figure 4(d), we show the close-up of the inner interface finger at the times 
$t=0.1118$, 
$0.1128$, 
$0.1138$, 
$0.1148$ and 
$0.1150$. At 
$t=0.1150$, the curvature at the cusp-like point is approximately 
$-274$, which is quite large compared with its initial value, 
$9.16\times 10^{-10}$. The distance between the inner interface and the sink is 
$1.19 \times 10^{-3}$, which is about twice the spatial resolution 
$\Delta x$. The GMRES iterations do not converge because of the resulting ill-conditioned linear system. As a note, we observe an excellent conservation of mass in the region 
$\varOmega _2(t)$. The area is preserved up to ten digits accuracy after the decimal point throughout the simulation.
In figures 5(a–c), because the distance between the two interfaces is small, the outer interface feels the presence of the sink quite strongly (though the sink is in 
$\varOmega _1$), and moves towards the sink along with the inner interface. By the time the inner interface starts to develop a finger to reach the sink, the outer interface is already very close to the inner interface, and the simulation stops at 
$T=0.1068$. In this simulation, at very early times, the distance between the two interfaces increases slightly from 0.0354 to 0.0357. That is because the inner interface moves a little faster than the outer one. Then the outer interface tends to catch up with the inner interface, and the distance between them decreases. The approaching velocity of the interfaces increases as time evolves, which follows approximately an exponential relation 
$0.0148\exp (43.39t)$. Unfortunately, we have to stop the calculation at later times, as the distance separating the two interfaces is quite small 
$2.05\times 10^{-3} \approx 3\,\Delta x$, where the spatial resolution 
$\Delta x$ is approximately 
$7.36\times 10^{-4}$. We found that the discretized linear system is very ill-conditioned, and the GMRES iteration solver does not converge.
These simulations suggest a new pattern forming mechanism by interface merging. The sink is the driving force; nevertheless, the precise nature of the instability is a result of the interaction between the two interfaces and the sink. As long as the interfaces are well separated, the inner interface approaches the sink before the outer one captures it; while if they are close initially, then it is more likely that they will come very close to each other before the inner interface reaches the sink.
4.5. Pattern formation with a sink and dissimilar outer and inner interfaces
The next extension is to consider cases where the inner and outer interfaces are no longer scaled versions of each other. We keep the outer interface as before (given by (4.5)), while the inner one is changed to a circle. All other parameters are the same as those used in figure 4. We set the radius of the inner interface, 
$r=0.7$ and 
$r=0.6$ initially, and display the simulation results in figures 6 and 7, respectively.
Effects of geometry. The sink is at the origin. The mobilities of the fluids are 
$M_1=1$, 
$M_2=100$ and 
$M_3=10\,000$. In (a), we display the initial configuration. At 
$t=0$, the outer interface is given by (4.5), and the inner interface is a circle with 
$r=0.7$. The intermediate and final configurations are shown in (b,c).
Effects of geometry. The sink is at the origin. The mobilities of the fluids are 
$M_1=1$, 
$M_2=100$ and 
$M_3=10\,000$. In (a), we display the initial configuration. At 
$t=0$, the outer interface is given by (4.5), and the inner interface is a circle with 
$r=0.6$. The intermediate and final configurations are shown in (b,c).
In figure 6(a), the outer and inner interfaces are placed quite close to each other on the negative 
$x$-axis at time 
$t=0$. At later times, the outer interface does not develop any fingers, while the inner interface shows an early sign of developing two fingers, marked by the two arrows pointing to the sink. However, before they fully develop, the outer interface comes very close to the inner one. The computation stops when the minimum distance between the interfaces is approximately 
$4.2 \times 10^{-4}$, while the spatial resolution 
$\Delta x$ is approximately 
$2.37 \times 10^{-4}$.
The development of these two fingers on the inner interface is far more prominent for 
$r=0.6$, as shown in figures 7(a–c). Here, we observe two well-developed fingers racing to the sink and forming two equal angles, giving the inner interface two distinct parts: (i) a bigger portion having a crescent shape; and (ii) a smaller region having the shape of an elongated drop along the negative 
$x$-axis. The inset of figure 7(c) zooms into the region where the parts of the inner interface come very close to each other. Compared with the single-interface case (Tian & Nie Reference Tian and Nie1998), the existence of a simple circular inner interface fundamentally alters the dynamics. At the end, the curvature at the cusp-like point is approximately 
$-250$. The distance between the inner interface and the sink is approximately 
$5.03 \times 10^{-3}$.
4.6. Pattern formation with a sink in the annulus region
We next consider that the sink is in the annular region, i.e. the fluid in 
$\varOmega _2$ gets extracted. We keep the outer interface the same as before. The inner interface is a circle of radius 
$r=0.2$ with its centre placed initially at (i) 
$(0.3,0)$, (ii) 
$(0.9,0)$ and (iii) 
$(-0.3,0)$. The sink remains at the origin. The results are summarized in figures 8(a–c), 9(a–c) and 10(a–c), where we show the evolution of morphologies of the two interfaces, and in figures 8(d), 9(d) and 10(d), where we plot the velocity of characteristic points on the interfaces at the positive (right) and negative (left) 
$x$-axis as a function of time.
Effects of sink location. The sink is at the origin. The mobilities of the fluids are 
$M_1=1$, 
$M_2=100$ and 
$M_3=10\,000$. In (a), we display the initial configuration. At 
$t=0$, the outer interface is given by (4.5), and the inner interface is a circle of radius 
$r=0.2$ with centre at 
$(0.3, 0)$. The intermediate and final configurations are shown in (b,c). In (d), we plot the velocity of various points on the interface.
Effects of sink location. The sink is at the origin. The mobilities of the fluids are 
$M_1=1$, 
$M_2=100$ and 
$M_3=10\,000$. In (a), we display the initial configuration. At 
$t=0$, the outer interface is given by (4.5), and the inner interface is a circle of radius 
$r=0.2$ with centre at 
$(0.9, 0)$. The intermediate and final configurations are shown in (b,c). In (d), we plot the velocity of various points on the interface.
Effects of sink location. The sink is at the origin. The mobilities of the fluids are 
$M_1=1$, 
$M_2=100$ and 
$M_3=10\,000$. In (a), we display the initial configuration. At 
$t=0$, the outer interface is given by (4.5). The inner interface is a circle of radius 
$r=0.2$ with centre at 
$(-0.3, 0)$. The intermediate and final configurations are shown in (b,c). In (d), we plot the velocity of various points on the interface.
In figures 8(c) and 9(c), we observe that the outer interface develops a finger with the tip on the negative 
$x$-axis, which moves towards the sink. It is evident that the inner interface is not circular, though in figure 9(c), it looks more circular. To quantify these results, we check the normal velocity of four characteristic points on the 
$x$-axis for both the inner and outer interfaces. As shown in figures 8(d) and 9(d), the point on the negative 
$x$-axis on the outer interface (in solid red) is the fastest moving. The near vertical segment of the curve implies the fact that the interface is rapidly approaching the sink towards the end of the simulation. On the other hand, the velocity of the point on the positive 
$x$-axis of the outer interface (in dashed red) is small, indicating that this point moves very slowly (normal velocity 
$\approx$ 0.01). The difference of the normal velocities inner left (in solid blue) and inner right (in dashed blue) in figure 8(d) explains the morphological change from the initial circular shape. In figure 9(d), we find that the normal velocities of both the points on the inner interface are nearly equal and very small, suggesting the better preservation of the circular shape of the inner interface.
For the case in figure 8(d), with our simulation data, the velocity seems to fit a relationship 
$0.7227(0.18499-t)^{-0.4654}+2.166$ even though in figure 9(d), the velocity seems to fit 
$4.126\times 10^{-4}(0.21671-t)^{-0.1.581}+4.639$. We note that our simulations stop at 
$t=0.1848$ and 
$t=0.21661$, respectively. Even though infinite velocity is not observed in our simulation, the velocity might blow up at a finite time. In figure 8(c), the distance between the outer interface and the sink is approximately 
$1.51 \times 10^{-3}$. For figure 9(c), the distance between the outer interface and the sink is approximately 
$2.63 \times 10^{-3}$, while the minimum distance between the interfaces is approximately 
$1.58 \times 10^{-2}$.
A significantly different scenario is observed when the centre of the inner domain is initially placed at 
$(-0.3, 0)$. As shown in figure 10(c), both interfaces are distorted considerably from their initial appearances. The region of the outer interface in the vicinity of the inner interface tends to bend around the inner one in its motion towards the sink, and eventually comes very close to the inner interface. This introduces a completely different final appearance for the outer interface as compared to the cases discussed above. The normal velocities also show a very different qualitative behaviour. The rightmost point on the outer interface moves quite fast towards the sink, while the leftmost point has a non-monotonic normal velocity. The velocity increases at early times when the outer and inner interfaces are well separated in space. But at later times, the velocity shows a rapid decrease when the outer interface comes close to the inner one. Note that the leftmost point on the inner interface hardly moves, while the rightmost point shows a normal velocity of magnitude close to 0.2. This novel ‘wrapping’ mechanism indicates non-trivial interactions between the two interfaces.
4.7. Effects of mobility
The relation between the viscosity ratio and the interfacial morphology is indeed a complicated one. In the case of a growing Hele-Shaw bubble, it is well known that at later stages of evolution, the fingering patterns depend strongly on the viscosity ratio of the fluids involved in a two-fluid system (Bischofberger, Ramachandran & Nagel Reference Bischofberger, Ramachandran and Nagel2015; Coutinho & Miranda Reference Coutinho and Miranda2020). Multi-layer cases have also been examined (Beeson-Jones & Woods Reference Beeson-Jones and Woods2015; Gin & Daripa Reference Gin and Daripa2015). However, the investigation of such problems in the case of sink-driven flow remains less explored. In this subsection, we check the effects of mobility in our problem.
In figures 11 and 12, we investigate the effects of mobility on the pattern formation. We select a new set of mobilities 
$M_1=0.01$, 
$M_2 =1$ and 
$M_3=100$ for regions 
$\varOmega _1$, 
$\varOmega _2$ and 
$\varOmega _3$, respectively, different from the values used before. We keep the capillary number 
$Ca$ unchanged. Therefore, other parameters remaining the same, the increase of mobilities by a factor of 100 can be understood as the same as increasing the rate of extraction 
$Q$ a hundred times, or decreasing the surface tension parameter by the same factor. In the initial configuration, the outer interface remains unchanged. We choose the inner interface to be a circle of radius 
$r=0.2$, centred at two different locations, 
$(-0.5, 0)$ and 
$(0.25,0)$. The sink stays at the origin and hence inside the annular region. In figure 11(c), which corresponds to the inner interface being placed at 
$(-0.5,0)$ initially, the outer interface tries to reach the sink from the left and almost wraps the inner interface. This results in the formation of two long fingers (fluid of 
$\varOmega _3$) penetrating into fluid in 
$\varOmega _2$ towards the sink, and a neck-like region of fluid 
$2$ wrapping the inner interface.
Effects of mobility. The sink is at the origin. The mobilities of the fluids are 
$M_1=0.01$, 
$M_2=1$ and 
$M_3=100$. In (a), we display the initial configuration. At 
$t=0$, the outer interface is given by (4.5). The inner interface is a circle of radius 0.2 initially centred at 
$(-0.5,0)$. The intermediate and final configurations are shown in (b,c). In (d), we show the tangent angle 
$\theta$ versus parametrization 
$\alpha$. In (e), we show the close-up of the outer interface finger at the times 
$t=0.1182$, 
$0.1198$ and 
$0.12133$.
Effects of mobility. The sink is at the origin. The mobilities of the fluids are 
$M_1=0.01$, 
$M_2=1$ and 
$M_3=100$. The outer interface is given by (4.5). The inner interface is a circle of radius 0.2 initially centred at 
$(0.25,0)$ in (a). In (b), we show the normal velocity of a left point on the outer interface, and in (c), the normal velocity of a left point on the inner interface, for various location of centres of the inner circle.
We suspect that at the tip of both fingers, two cusp-like singularities are about to form. This is because, in the adjacent 
$\theta \unicode{x2013} \alpha$ plot of the outer interface shown in figure 11(d), where 
$\theta$ is the tangent angle and 
$\alpha$ is the parameter that parametrizes the outer interface, we observe a sharp transition in the value of 
$\theta$ near the region marked ‘a’ in the inset of figure 11(c), similar to the calculation of cusp-like formation observed earlier (Tian & Nie Reference Tian and Nie1998). Figure 11(e) shows a close-up of the inner interface fingers at the times 
$t=0.1182$, 
$0.1198$ and 
$0.12133$.
A remarkable situation is observed when the centre of the inner interface is placed at 
$(0.25,0)$ at 
$t=0$, shown in figure 12(a). Here, we plot a sequence of morphologies as insets, and the curves are the 
$\theta \unicode{x2013} \alpha$ relation of the inner interface at different times close to the point where simulation fails. The inner interface clearly experiences the presence of the sink and is pulled strongly towards it, forming what looks to be a small but distinct cusp-like pattern, which distinguishes itself from those reported earlier by Tian & Nie (Reference Tian and Nie1998): the cusp-like morphology in our case forms in the outward direction, whereas in their case the cusp is inwards. Again, the 
$\theta \unicode{x2013} \alpha$ curve of the inner interface shows a steep transition near 
$\alpha = 3.14$. Finally, the normal velocity plots of the left point on the outer and inner interfaces are shown in figures 12(b) and 12(c), respectively, for various initial locations of the initial centre of the inner circle. The curves indicate that the motion of the two interfaces approaches a possible (quasi-)steady state as time progresses.
4.8. Evolution of surface energy
In this subsection, we investigate the evolution of surface energy. The energy of the interface is defined as 
$E(t)= \int _{\varGamma (t)} \sigma \,\textrm {d}s$, where 
$\sigma$ is the surface tension of the interface 
$\varGamma (t)$. Apparently, 
$E(t)$ is related to the length of the interface. For example, when the sink is at the centre of a circle, 
$E(t)=2\sigma \sqrt {{\rm \pi} \,A(t)}$ in theory for a constant surface tension, where 
$A(t)$ is the area enclosed by the interface at time 
$t$. To scale out the size of the interface, we consider a non-dimensional energy and obtain 
$E(t)/E(0)=\sqrt {A(t)/A(0)}$ for the circle. When all the fluid is removed, the energy goes to zero, i.e. the interface shrinks to a point.
First, we explore the evolution of 
$E(t)/E(0)$ with respect to 
$A(t)/A(0)$ for the inner interface in our simulations when the sink is in the interior region of the inner interface. Different types of various initial configurations of the two interfaces are used (see table 1). The result is summarized in figure 13(a), and one dataset is chosen from each type. The theoretical formula is used for the result of two circles, where the sink is placed right at the centre of the annulus. Our simulation results coincide with the formula, while the area does not go to zero in our simulation. For other cases, the interface remains smooth and the energy decreases as the area shrinks at early times. When the inner interface experiences long fingers, the energy starts to increase while the area decays. The energy evolution of simulations shown in § 4.4 (figure 4) and § 4.5 (figure 7) also demonstrate the behaviour. Note that for the case in figure 7, the inner interface is a circle centred at the origin. The interactions between the inner interface and outer interface make the inner interface deform from the circle. It takes a long time for the interface to form long fingers. Thus the energy evolution agrees very well with the theoretical results (two circles) for a long period.
Initial configurations of the two interfaces when the sink is inside the inner interface.
(a) The evolution of non-dimensional energy 
$E(t)/E(0)$ with respect to the non-dimensional area 
$A(t)/A(0)$ when the sink is in the inner interface. (b) The evolution of non-dimensional energy 
$E(t)/E(0)$ with respect to the non-dimensional area 
$A(t)/A(0)$ when the sink is in the annulus region.
Next, we investigate the evolution of 
$E(t)/E(0)$ for the outer interface when the sink is in the fluid region 
$\varOmega _2$, i.e. 
$A(t)$ represents the area of fluid domain 
$\varOmega _2$, and 
$E(t)$ is the surface energy of 
$\varGamma _2(t)$. Again, different types of various initial configurations of the two interfaces are used (see table 2). The result is summarized in figure 13(b), and one dataset is chosen from each type. At early times, the energy and the area both decay. These data could be fitted as 
$E(t)/E(0)\sim (A(t)/A(0))^{0.658}$. During this moment, the outer interface is smooth and does not experience long fingers. Later, multiple long fingers are formed and the surface energy starts to grow. The energy evolution of simulations shown in § 4.6 (figures 8 and 10) and § 4.7 (figure 11) also satisfies the behaviour.
Initial configurations of the two interfaces when the sink is in the annulus region.
5. Conclusion
In this paper, we have investigated a three-layer Hele-Shaw problem where the interfaces move due to the presence of a sink. We present the governing equations of the problem and the corresponding boundary integral formulation. The boundary integral equations are discretized by spectrally accurate quadratures, and we march in time with a second-order-accurate time-stepping technique after alleviating the numerical stiffness issue. We have performed simulations by varying the initial geometry of the interfaces, the location of sink and the mobility of the fluids. In a single-interface problem, the singularity occurs when the interface reaches the sink. However, in the multi-layer problem, we observe that the singularity may also occur because the interfaces come very close to each other. We observe rich interface dynamics, and report novel cases beyond those reported previously in the literature (Tian & Nie Reference Tian and Nie1998). A natural extension of our work would be to consider a more practical geometry consisting of multiple inner regions of fluids 
$\varOmega _{11}, \varOmega _{12}, \ldots, \varOmega _{1n}$ instead of just 
$\varOmega _{1}$, which would be all surrounded by the interface 
$\varGamma _2$. This would better capture the scenario where multiple air bubbles (or regions of less viscous fluids) are trapped within, for example, a system with two additional fluids. Further, it would be interesting to see the effects of multiple sinks and, possibly, a combination of both sources and sinks.
Apart from the Hele-Shaw community, the current work could be of interest for the multi-phase flows in permeable media, where the flow is well approximated by Darcy's laws. In oil-filled rocks, several fluids such as water, oil and air might be present. There could be areas where all these fluids of contrasting viscosity flow while interacting with each other. Our present work would prove quite relevant in such situations, especially if we extend this problem to more complicated set-ups.
Finally, we make some comments on the experimental side of our problem. It is not hard to find experimental studies looking into a single-layer problem (Logvinov Reference Logvinov2019). Multi-layer experiments are few (Cardoso & Woods Reference Cardoso and Woods1995); however, these studies are gaining momentum. For example, there has been work to explore the effect of squeezing the plates, which reduces the gap distance, on the flow pattern (Moffatt, Guest & Huppert Reference Moffatt, Guest and Huppert2021). A very recent work explores the gravity-driven flow in four-layer cells (Brahim & Thoroddsen Reference Brahim and Thoroddsen2022). Since the external force, i.e. gravity, drives the flow, the arrangement does not need injection or removal of the fluid. Finally, we wish to add that suction-related Hele-Shaw experiments are limited, and our findings might serve as a benchmark for the experimental fluid mechanics community, and shed light on or motivate the development of analytical works for this practically important problem.
Funding
S.L. and J.S.L. acknowledge support from the National Science Foundation, grants DMS-1720420, DMS-2309800 and DMS-2309798. W.Y. and M.Z. acknowledge support from the National Science Foundation of China, grants DMS-11771290 and DMS-12301553. M.Z. also thanks for support the Ministry of Education Key Lab in Scientific and Engineering Computing.
Declaration of interests
The authors report no conflict of interest.
Appendix
In this Appendix, we display the initial configurations corresponding to figures 13(a) and 13(b) in tables 1 and 2, respectively.
