1. Introduction
Viscous fingering (VF) or Saffman–Taylor instability (Saffman & Taylor Reference Saffman and Taylor1958) occurs when a less viscous fluid displaces a more viscous fluid in porous media, Hele-Shaw cells or microfluidic cells. Since it was first reported by Hill (Reference Hill1952), extensive research investigations have been performed due to its widespread applications (Homsy Reference Homsy1987; De Wit Reference De Wit2020). Most of the reported works, however, are mainly on a single displacement front. Interests given to dual displacement fronts, in which two fluid interfaces are formed, are still limited despite their importance for various applications. Given the characteristics of dual displacement fronts, it can be harnessed in two opposite ways. The first is to generate efficient fluid displacement processes in reservoir treatment (Paraskeva et al. Reference Paraskeva, Charalambous, Stokka, Klepetsanis, Koutsoukos, Read, Ostvold and Payatakes2000; Talaghat, Esmaeilzadeh & Mowla Reference Talaghat, Esmaeilzadeh and Mowla2009), enhanced oil recovery (Le Van & Chon Reference Le Van and Chon2017; Vishnudas & Chaudhuri Reference Vishnudas and Chaudhuri2017; Afzali, Rezaei & Zendehboudi Reference Afzali, Rezaei and Zendehboudi2018; Chaudhuri & Vishnudas Reference Chaudhuri and Vishnudas2018; Bakharev et al. Reference Bakharev, Enin, Groman, Kalyuzhnyuk, Matveenko, Petrova, Starkov and Tikhomirov2022), in situ groundwater remediation (Wood, Simmons & Hutson Reference Wood, Simmons and Hutson2004) and column chromatography (Mayfield et al. Reference Mayfield, Shalliker, Catchpoole, Sweeney, Wong and Guiochon2005; Shalliker et al. Reference Shalliker, Catchpoole, Dennis and Guiochon2007; Shalliker & Guiochon Reference Shalliker and Guiochon2010) for hydrogeology, petroleum and chemical processes. The second is to induce fluid mixing in a confined fluid system with a small Reynolds number for microdroplet generation (Cardoso & Woods Reference Cardoso and Woods1995; Hashimoto et al. Reference Hashimoto, Garstecki, Stone and Whitesides2008) and microflow mixing (Glasgow & Aubry Reference Glasgow and Aubry2003; Coleman & Sinton Reference Coleman and Sinton2005; MacInnes, Chen & Allen Reference MacInnes, Chen and Allen2005) for various biological and chemical syntheses. To generate efficient fluid displacement, stable displacement should be maintained by avoiding VF (Yuan & Azaiez Reference Yuan and Azaiez2014; Sharma et al. Reference Sharma, Bin Othman, Nagatsu and Mishra2021). In contrast, VF is desirable for enhancing the mixing process by breaking the stability (Jha, Cueto-Felgueroso & Juanes Reference Jha, Cueto-Felgueroso and Juanes2013; Chen et al. Reference Chen, Huang, Huang and Miranda2015). To the best of our knowledge, the exploration of these strategies is still an open question.
The first work on the dual displacement fronts was performed by De Wit, Bertho & Martin (Reference De Wit, Bertho and Martin2005), and afterwards, a series of similar studies were reported (Mishra, Martin & De Wit Reference Mishra, Martin and De Wit2008, Reference Mishra, Martin and De Wit2009, Reference Mishra, Martin and De Wit2010; Hota, Pramanik & Mishra Reference Hota, Pramanik and Mishra2015). Some experimental studies have also been performed by utilizing micropillar array columns (De Malsche et al. Reference De Malsche, Op De Beeck, Gardeniers and Desmet2009; Haudin et al. Reference Haudin, Callewaert, De Malsche and De wit2016). However, the analyses were mainly focused on the VF structure, the linear stability analysis and the mixing process, whereas analysis of the interaction dynamics of VF has never been discussed. In this paper, we explored the interaction dynamics of VF in dual displacement fronts through numerical simulation by varying the viscosity configuration. By using each interface onset time, finger velocity and the finite slice breakthrough time, we modelled and mapped the interaction dynamics of VF in dual displacement fronts with various viscosity configurations of the finite slice and bulk fluid.
2. Mathematical formulation and numerical solution
We consider miscible fluids 
$A$, 
$B$ and 
$C$ placed from left to right, forming two-fluid interfaces of 
$AB$ and 
$BC$, referred to as the rear interface and the front interface, respectively. The fluid flow is considered in a two-dimensional porous medium with uniform porosity and constant permeability 
$\kappa$. Fluid 
$A$ then flows to the right, pushing the dual displacement front configuration, as shown in figure 1.
Figure 1. Sketch of the fluid configuration system.
With the assumptions that the system is incompressible, neutrally buoyant and first contact miscible, the governing equations are
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\underline{u} = 0, \quad \boldsymbol{\nabla}p ={-} \displaystyle\frac{\mu}{\kappa}\underline{u}, \end{gather}$$
$$\begin{gather}\displaystyle\frac{\partial a}{\partial t} + \underline{u}\boldsymbol{\cdot}\boldsymbol{\nabla}a = D_{A}\nabla^2a, \end{gather}$$
$$\begin{gather}\displaystyle\frac{\partial b}{\partial t} + \underline{u}\boldsymbol{\cdot}\boldsymbol{\nabla}b = D_{B}\nabla^2b, \end{gather}$$
$$\begin{gather}\displaystyle\frac{\partial c}{\partial t} + \underline{u}\boldsymbol{\cdot}\boldsymbol{\nabla}c = D_{C}\nabla^2c. \end{gather}$$ In the above expressions, the fluid flow velocity 
$\underline {u}=(u,v)$, 
$p$ is the pressure, 
$\mu$ is the fluid viscosity, 
$t$ is time, 
$a$, 
$b$ and 
$c$ are the concentration of the fluids 
$A$, 
$B$ and 
$C$, respectively. It is assumed that the diffusion coefficients of fluids 
$A$, 
$B$ and 
$C$ are the same, i.e. 
$D_A=D_B=D_C=D$. The initial concentration value of 
$a$ is 
$a_{0}$, and subscripts 
$A$, 
$B$ and 
$C$ represent the parameters of fluids 
$A$, 
$B$ and 
$C$, respectively. The viscosity concentration relation is considered as Arrhenius type (Hejazi et al. Reference Hejazi, Trevelyan, Azaiez and De Wit2010; Hejazi & Azaiez Reference Hejazi and Azaiez2012),
\begin{equation} \mu = \mu_A \exp\left({\left(R_{B}\displaystyle\frac{b}{a_0} + R_{C}\frac{c}{a_0}\right)}\right), \end{equation}
where, 
$R_{B} = \ln (\mu _{B}/\mu _{A})$, and 
$R_{C} = \ln ({\mu _{C}}/{\mu _{A}})$ are the log mobility ratios between the viscosity of 
$B$, 
$C$ with 
$A$, respectively. The governing equations above were then non-dimensionalized by using injection velocity (
$U$), hydrodynamical time (
$\tau _h = {D}/{U^2}$) and length (
$L_h = {D}/{U}$). Viscosity and pressure are normalized by 
$\mu _A$ and 
$\mu _A D/ \kappa$, respectively. The concentration of species 
$A$, 
$B$ and 
$C$ are also scaled by 
$a_0$ (initial concentration of 
$A$). A reference frame moving with the injection velocity was used to focus on the dual displacement front dynamics. Afterwards, introducing a stream function 
$\psi (x,y)$ such that 
$u = {\partial \psi }/{\partial y}$ and 
$v = -{\partial \psi }/{\partial x}$, the model equations can be written in the stream function–vorticity formulation (Tan & Homsy Reference Tan and Homsy1988; De Wit et al. Reference De Wit, Bertho and Martin2005), as
$$\begin{gather} \nabla^2{\psi} ={-} R_B(\psi_x b_x + \psi_y b_ y + b_y) - R_C(\psi_x c_x + \psi_y c_ y + c_y), \end{gather}$$
$$\begin{gather}a_t + a_x \psi_y - a_y \psi_x = \nabla^2a, \end{gather}$$
$$\begin{gather}b_t + b_x \psi_y - b_y \psi_x = \nabla^2b, \end{gather}$$
$$\begin{gather}c_t + c_x \psi_y - c_y \psi_x = \nabla^2c, \end{gather}$$
$$\begin{gather}\mu=e^{R_B b + R_C c}. \end{gather}$$ For the numerical simulation, dimensions of 
$L_y=1024$, 
$L_x=5120$ and 
$W=1024$ (figure 1) were used. This slice width was selected because it gives enough space for the fingering development without the disturbance from the periodic boundary layer. In addition, it gives a reasonable computational time to observe the breakthrough time. We consider periodic boundary conditions in both the axial and transverse directions to avoid disturbances to the fluid flow. They do not affect the dual displacement front dynamics at the centre. An initial random disturbance was introduced to induce the instability. As reported by De Wit et al. (Reference De Wit, Bertho and Martin2005), this random disturbance magnitude mainly affects the VF onset time. Smaller random disturbance magnitude leads to longer onset time, resulting in unnecessarily longer computational time. Therefore, to ensure similar initial conditions and reasonable computational time, a constant value of random disturbance of 
$O(10^{-2})$ was selected throughout this work. Further, (2.6)–(2.10) are numerically solved using a Fourier pseudospectral method introduced by Tan & Homsy (Reference Tan and Homsy1988), which has been shown to successfully model various numerical studies of VF (De Wit & Homsy Reference De Wit and Homsy1999; Nagatsu & De Wit Reference Nagatsu and De Wit2011).
For the quantitative evaluation, parameters based on the transverse average fluid concentration profile
\begin{equation} { i_{avg,y}(x,t) = \frac{1}{L_y} \int_{0}^{L_y}i(x,y,t)\,{{\rm d}y} }, \end{equation}
with 
$i$ corresponding to the concentrations (
$a$, 
$b$ and 
$c$) of fluid 
$A$, 
$B$ or 
$C$ were used (see figure 2). First, the mixing length (
$ML$) in red of the fluids 
$A$ and 
$C$ defined as the length between the transverse average fluid concentration of 
$\epsilon$ and 
$(\epsilon -1)$ with 
$\epsilon =0.01$ was obtained. Next, the deformation of the rear interface (
$\tau _{ons,R}$) and the front interface (
$\tau _{ons,F}$) were measured by comparing the mixing length with the pure diffusion case. The onset time is determined as the time when the mixing length starts deviating from the pure diffusive interface value by 10. Lastly, we introduced the breakthrough time 
$\tau _{bk}$ that was determined when any transverse average fluid concentration contained concentrations of both fluid 
$A$ and 
$C$ larger than 
$\epsilon$, indicating the breakthrough of fluid 
$B$ slice (see figure 2). Once the fluid 
$B$ slice has broken through, the dual displacement fronts become inefficient for displacement (Paraskeva et al. Reference Paraskeva, Charalambous, Stokka, Klepetsanis, Koutsoukos, Read, Ostvold and Payatakes2000; Wood et al. Reference Wood, Simmons and Hutson2004; Mayfield et al. Reference Mayfield, Shalliker, Catchpoole, Sweeney, Wong and Guiochon2005; Shalliker et al. Reference Shalliker, Catchpoole, Dennis and Guiochon2007; Talaghat et al. Reference Talaghat, Esmaeilzadeh and Mowla2009; Le Van & Chon Reference Le Van and Chon2017; Vishnudas & Chaudhuri Reference Vishnudas and Chaudhuri2017; Afzali et al. Reference Afzali, Rezaei and Zendehboudi2018; Chaudhuri & Vishnudas Reference Chaudhuri and Vishnudas2018; Bakharev et al. Reference Bakharev, Enin, Groman, Kalyuzhnyuk, Matveenko, Petrova, Starkov and Tikhomirov2022) but effective for fluid mixing (Coleman & Sinton Reference Coleman and Sinton2005; MacInnes et al. Reference MacInnes, Chen and Allen2005; Jha et al. Reference Jha, Cueto-Felgueroso and Juanes2013; Chen et al. Reference Chen, Huang, Huang and Miranda2015). Therefore, this parameter is used as the main indicator in characterising the dual displacement fronts instability dynamics.
Figure 2. Sketch of the onset and breakthrough concepts.
3. Results and discussion
Simulations with 
$R_C=2$, 
$R_C=0$ and 
$R_C=-1$, representing the condition of 
$R_C>0$, 
$R_C=0$ and 
$R_C<0$ were performed with varied 
$R_B$. This helps to compare the results when the viscosity of displacing fluid 
$A$ is more, the same or less than the viscosity of displaced fluid 
$C$. Simulations are performed for five different random numbers generated for the seeding of the instability in each set of parameters. In order to find the onset of deformation, we computed the evolution of mixing length for several 
$R_B$ and 
$R_C=2$, 0 and 
$-1$, with the same initial condition and plotted in figure 3. The examples of mixing length development for all 
$R_B$ and 
$R_C$ at one randomness case are shown. The onset of instabilities of both frontal–front and rear–front can be seen in figures 3(a–c) and 3(d–f), respectively, when the mixing length departs from the respective pure diffusive mixing lengths (stable case). As seen from figure 3(a), for 
$R_C=2$, and for increasing of 
$R_B$ from 
$-1$ until 
$R_B=0.7$ onset of instability at the 
$AB$ front delayed and after a critical 
$R_B \gtrsim 0.8$ (star-marked curve in figure 3a), the onset of instability reversed and became early. A similar trend occurs at the 
$BC$ front as depicted in figure 3(d), i.e. the reversal of delayed onset to early onset occurs when 
$R_B \gtrsim 1.4$ (star-marked curve in figure 3d). However, for 
$R_C=0$ and 
$-1$ when log mobility ratios 
$|R_B|$ increase and get farther away from 
$R_A$ and 
$R_C$, the early onset is observed, and no reversal behaviour is seen on the onset of instability.
Figure 3. Mixing length of fluid (a–c) 
$A$ and (d–f) 
$C$ at 
$R_C$ of 2, 0 and 
$-1$ for various values of 
$R_B$ from one of the randomness case. The bottom black dashed line shows the pure diffusive mixing length of stable interfaces.
Further, the numerical simulation results for 
$R_C=2$ with 
$R_B$ of 
$-0.4$, 0.4, 1.0, 1.6 and 2.4 are depicted in figure 4 through the density plot of concentrations of species B. Initially, the frontal interface was unstable (e.g. see Supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.670, with 
$R_B=0.4, R_C=2.0$), followed by the dual front (e.g. Supplementary movie 2, with 
$R_B=1.0, R_C=2.0$) and then the rear interface became unstable (e.g. Supplementary movie 3, with 
$R_B=1.6, R_C=2.0$). At a specific time, the plots are shown representing the onset of VF front 
$\tau _{ons,R}$, 
$\tau _{ons,F}$ and breakthrough time 
$\tau _{bk}$, which were chosen from the mixing length shown in figure 3.
Figure 4. Numerical simulation results of 
$R_C=2$ with various value 
$R_B$ corresponding to the (a–e) 
$\tau _{on,R}$ or 
$\tau _{on,F}$ and 
$\tau _{bk}$ (f–j). Here 
$R_B$ of 
$-0.4$ and 0.4 demonstrated a front-dominated fingering, whereas 
$R_B$ of 1.6 and 2.4 demonstrated a rear-dominated fingering. Here 
$R_B$ of 1.0, on the other hand, demonstrated a dual fingering.
With 
$R_B$ of 
$-0.4$ (figure 4a,f), the rear interface is stable, whereas the front interface is unstable. As a result, the VF developed at the front interface and then disturbed the stable rear interface. For the 
$R_B$ of 0.4 (figure 4b,g), similar phenomena with 
$R_B$ of 
$-0.4$ also occur, but later. Even though both interfaces are unstable, the rear interface is more stable than the front interface due to the higher viscosity contrast between fluids 
$B$ and 
$C$ than fluids 
$A$ and 
$B$. As a result, VF at the front interface occurred earlier and disturbed the rear interface before the instability grew. Similar but opposite phenomena also occur with 
$R_B$ of 1.6 and 2.4 (figure 4d,e,i,j), in which the VF at the rear interface grew earlier and disturbed the front interface. These four cases demonstrate a single instability case with the 
$R_B$ of 
$-0.4$ and 0.4 exhibiting front-dominated fingering and 
$R_B$ of 1.6 and 2.4 exhibiting rear-dominated fingering, respectively. Although in 
$R_B$ of 0.4 and 1.6, both interfaces are unstable, the more unstable interface disturbs the less stable interface before the instability at the less stable interface is developed. On the other hand, in the case 
$R_B$ of 1.0 (figure 4c,h), both interfaces are equally unstable from the double instability. As a result, the VF started growing at the same time and finally collided, resulting in dual fingering. These instability dynamics can be divided into two types. The first type is the ordinary VF initiated by the interface viscosity contrast. The second type is not an instability (because the interface is stable) but rather a deformation initiated from VF of the other interface. Therefore, for this type, the onset time and the finger velocity of the other interface govern the onset of this instability mechanism.
Plotting 
$\tau _{ons,R}$, 
$\tau _{ons,F}$ and 
$\tau _{bk}$ with 
$R_B$ in figure 5 provides a comprehensive point of view on the instability dynamics. By gradually increasing 
$R_B$ from 
$-1.0$, 
$\tau _{ons,F}$ increases because the front interface is becoming more stable. Because this condition demonstrates front-dominated fingering, 
$\tau _{ons,R}$ and 
$\tau _{bk}$ coincide, in which they also increase with 
$R_B$ due to the delay in 
$\tau _{ons,F}$. When 
$R_B$ is higher than 0.7, the delay on 
$\tau _{ons,F}$ provides sufficient time for the rear interface to develop instability. As a result, 
$\tau _{ons,R}$ becomes smaller than 
$\tau _{bk}$, commencing the dual fingering regime. As 
$\tau _{ons,F}$ is getting delayed, when 
$R_B$ is 1.0, 
$\tau _{ons,R}$ and 
$\tau _{ons,F}$ occur at relatively the same time. However, at 
$R_B$ of 1.3, 
$\tau _{ons,F}$ and 
$\tau _{bk}$ start coinciding, marking the beginning of the rear-dominated fingering regime. Past this point, 
$\tau _{ons,R}$ keeps decreasing because the rear interface is getting less stable, and thus 
$\tau _{ons,F}$ and 
$\tau _{bk}$ also keep decreasing. Given these fingering dynamics, the 
$\tau _{bk}$ of 0.7 and 1.3 can be identified as the critical fluid 
$B$ viscosity 
$R_{B,crit,1}$ and 
$R_{B,crit,2}$, in which the regime transition occurs.
Figure 5. The 
$\tau _{ons,R}$, 
$\tau _{ons,F}$ and 
$\tau _{bk}$ of 
$R_C=2$ at various 
$R_B$. Data points correspond to the averaged value from five simulations, and the colour-shaded area represents the data scattering of the five simulations. The solid lines represent the model from (3.1), (3.2) and (3.6). The green and the cyan dashed lines correspond to 
$\mu _A=\mu _B$ and 
$\mu _B=\mu _C$, respectively, whereas the red dashed lines correspond to both 
$R_{B,crit,1}$ and 
$R_{B,crit,2}$.
These instability dynamics maps can be modelled by first generating the model for the onset time of both interface (
$\tau _{ons,R}$ and 
$\tau _{ons,F}$), diffusion propagation (
$X_D$) and both finger downstream and upstream velocities (
$V_+$ and 
$V_-$). For the onset time, it decreases with the interface viscosity ratio by following the power function of minus two (Tan & Homsy Reference Tan and Homsy1988; De Wit et al. Reference De Wit, Bertho and Martin2005) (figure 6a) as follows:
$$\begin{gather} \tau_{ons,R} = 2215R_B^{{-}2} \quad {\rm for} \ R_B> 0, \end{gather}$$
$$\begin{gather}\tau_{ons,F} = 2215(R_C-R_B)^{{-}2} \quad {\rm for}\ R_C-R_B> 0. \end{gather}$$
These models fit well with the data with a coefficient of determination (
$r^2$) of 0.991 (see figure 6a). The green and blue solid lines in figure 5 clearly depict these models. The propagation of the diffusion fronts is modelled with a time square-root model by using the data of interface propagation at a stable state (figure 6b). The data fit well with a coefficient of determination (
$r^2$) of 0.993. The corresponding model for the case of a pure diffusion front is given in (3.3) with the proportionality constant 
$\alpha$ as 3.29 when 
$\epsilon$ is chosen as 0.01. For the finger velocity, the upstream and downstream directions need to be modelled separately because the upstream finger velocity tends to be slower than the downstream finger velocity (Mishra et al. Reference Mishra, Martin and De Wit2010; Nijjer, Hewitt & Neufeld Reference Nijjer, Hewitt and Neufeld2018) (figure 6c). The finger velocity was calculated by measuring the required time for the finger to reach the other interface from its onset to its breakthrough, which can be obtained in the regimes of the front-dominated fingering and rear-dominated fingering by taking into account the diffusion propagation as well. Because the diffusion propagation is independent of viscosity ratio but the finger velocity is a function of viscosity ratio (Nijjer et al. Reference Nijjer, Hewitt and Neufeld2018), the following models were generated:
$$\begin{gather} X_D(t) = X_\epsilon = X_{1-\epsilon} = \alpha t^{0.5} \quad {\rm with}\ \alpha = f(\epsilon) =3.29, \end{gather}$$
$$\begin{gather}V_+= 0.285R_B \quad {\rm for} \ R_B> 0, \end{gather}$$
$$\begin{gather}V_-= 0.197(R_C-R_B) \quad {\rm for}\ R_B> 0, \end{gather}$$
where 
$X_\epsilon$ and 
$X_{1-\epsilon }$ are the locations of the maximum and minimum threshold concentration for mixing length, which is at 
$\epsilon = 0.01$. From these model equations ((3.3)–(3.5)), the breakthrough time 
$\tau _{bk}$, as depicted by the solid red line in figure 5 on the three regimes of fingering dynamics, can be modelled with an equation as follows:
\begin{equation} \tau_{bk} = \frac{W-X_D(t) +\tau_{ons,R}V_++\tau_{ons,F}V_-}{V_++V_-}\end{equation}with
$$\begin{gather} V_+= 0 \quad {\rm if} \ \tau_{ons,R} \geqslant \frac{W-X_D}{V_-}+\tau_{ons,F} \quad {\rm at}\ R_{B,crit,1}, \end{gather}$$
$$\begin{gather}V_-= 0 \quad {\rm if} \ \tau_{ons,F} \geqslant \frac{W-X_D}{V_+}+\tau_{ons,R} \quad {\rm at}\ R_{B,crit,2}, \end{gather}$$
$$\begin{gather}\tau_{bk} = \left(\frac{W}{2\alpha}\right)^2 \quad {\rm if} \ W \leqslant 2X_D. \end{gather}$$
Figure 6. (a) The model formulation for 
$\tau _{ons,R}$ and 
$\tau _{ons,F}$ as shown in (3.1) and (3.2). (b) The model formulation for 
$X_D$ as in (3.3). (c) The model formulation for (
$V_+$) and (
$V_-$) as shown in (3.4) and (3.5) and (d,e) the iteration process for 
$\tau _{bk}$ at (3.6)–(3.9) for 
$R_C$ of 2 and 0.
Given the time variable presented in 
$X_D$ at (3.3), which also contributes to breakthrough time 
$\tau _{bk}$, hence, we calculate the breakthrough time 
$\tau _{bk}$ given in (3.6) using an iterative process with 
$X_D(t)$ given in (3.3) and defined as
$$\begin{gather} \tau_{bk}^{j} = \displaystyle\frac{W-\alpha (t^j)^{0.5} +\tau_{ons,R}V_++\tau_{ons,F}V_-}{V_++V_-}, \end{gather}$$
$$\begin{gather}\tau_{bk}^{j+1} = \displaystyle\frac{W-\alpha ( \tau_{bk}^{j} )^{0.5} +\tau_{ons,R}V_++\tau_{ons,F}V_-}{V_++V_-}. \end{gather}$$ We performed iterations considering the initial time as 
$t=0$, and with three iterations, we achieved the convergence with a relative error of 
$O(10^{-4})$. The result of the breakthrough time obtained from the iteration process is shown in figures 6(d) and 6(e) for 
$R_C$ of 2 and 0. We found that the value of 
$\tau _{bk}$ only changes slightly during each iteration. As for the pure diffusion condition given in (3.9), since it is only a function of the diffusion coefficient with a constant slice width of 1024, the 
$\tau _{bk}$ value is always fixed at 24 219 as also shown in figure 6(e).
As shown in figure 5, these models can predict the dynamics accurately. The 
$\tau _{ons,R}$ follows the model for 
$R_B \geqslant R_{B,crit,1}$, and the 
$\tau _{ons,F}$ follows the model for 
$R_B \leqslant R_{B,crit,2}$. The onset time beyond this condition no longer follows the model because the instability is not governed by viscosity contrast but by the disturbance from the other unstable interface. These models can also predict both 
$R_{B,crit,1}$ and 
$R_{B,crit,2}$ as the boundary between the three fingering regimes. The higher peak of 
$\tau _{bk}$ located at 
$R_{B,crit,1}$ corresponds to the slower 
$V_-$ than the 
$V_+$, leading to later disturbance to the rear interface. The least accurate prediction is in the dual-fingering regime. Because the finger velocity models are based on the fully developed finger, these models cannot catch the early nonlinear finger dynamics near the onset. Given that two fingers were developing, this effect becomes more significant.
Now, for 
$R_C=0$, the stability can also be divided into three regimes (figure 7a), but instead of dual fingering regime, a stable regime was found. Supplementary material (movie 4) provides the animations of stable flow dynamics with the condition of 
$R_B = -0.5$ and 
$R_C = -1.0$, which also helps to compare the different unstable regimes. At 
$R_B<0$, the rear interface is stable, but the front interface is unstable, resulting in front-dominated fingering, whereas 
$R_B>0$ shows the opposite, in which the rear interface is unstable, but the front interface is stable, corresponding to the rear dominated fingering. At 
$R_B=0$, on the other hand, both interfaces are stable, resulting in a stable regime without 
$\tau _{ons,R}$ and 
$\tau _{ons,F}$. However, 
$\tau _{bk}$ can still occur due to the diffusion propagation, although at a much slower rate. Therefore, with 
$R_B$ closing to 0, 
$\tau _{bk}$ increase sharply until it reaches a plateau of pure diffusion, spanning beyond the stable regime because breakthrough by diffusion propagation is still possible when the 
$\tau _{ons,R}$ and 
$\tau _{ons,F}$ occur late. This also serves as the 
$\tau _{bk}$ limit of the dual displacement system when 
$\tau _{bk}$ occurs due to diffusion propagation instead of VF interaction, as also given in (3.6). For 
$R_C=-1$, instability dynamics similar to those of 
$R_C=0$ were also observed (figure 7b). The generated models of 
$\tau _{ons,R}$, 
$\tau _{ons,F}$ and 
$\tau _{bk}$ also demonstrate satisfying agreement with the data. The only difference is that the stable regime expands to 
$R_B$ between 0 and 
$-1$.
Figure 7. Onset and breakthrough time, 
$\tau _{ons,R}$, 
$\tau _{ons,F}$ and 
$\tau _{bk}$ of (a) 
$R_C = 0$ and (b) 
$R_C= - 1$ at various 
$R_B$. Data points correspond to the average value of five simulations with the standard deviation. The solid lines represent the generated model from (3.1), (3.2) and (3.6). The green and the cyan dashed line correspond to 
$\mu _A=\mu _B$ and 
$\mu _B=\mu _C$, respectively, whereas the purple area corresponds to the stable regime.
Based on the developed models, the regime map complemented with the prediction of 
$\tau _{bk}$ can be created by extending the calculation for various 
$R_B$ and 
$R_C$ as shown in figure 8. Four regimes of rear-dominated fingering, dual fingering, front dominated fingering and stable are given. A stable regime is bound by the lines of 
$R_B=R_C$ and 
$R_B=0$ and surrounded by the region of pure diffusion, as shown by the white-coloured area. The dual fingering regime is bound by the lines of 
$R_{B,crit,1}$ and 
$R_{B,crit,2}$ as given in (3.7) and (3.8). The remaining area corresponds to the front-dominated fingering at the left and rear-dominated fingering at the right.
Figure 8. The regime map complemented with the 
$\tau _{bk}$ prediction for 
$R_B$ of 
$-5$ to 5 and 
$R_C$ between 
$-5$ and 5. The black dashed lines correspond to the regime boundary. The pink dotted lines correspond to the 
$R_C=2$, 
$R_C=0$ and 
$R_C=-1$, whereas the brown dotted lines correspond to 
$R_B=-1$ and 
$R_B=0.8$.
In addition, further classification can be performed by fixing either 
$R_B$ or 
$R_C$. The system with 
$R_C<0$ is inherently stable because the system can become stable by changing 
$R_B$ in between 
$R_C$ and 0. In contrast, for 
$R_C>0$, the configuration is inherently unstable because the system can never be stable unless the viscosity differences are very small, resulting in more dominant diffusion propagation than fingering. Similarly, by changing 
$R_C$, 
$R_B<0$ is also inherently stable because it can be stable when 
$R_C$ is lower than 
$R_B$, whereas 
$R_B>0$ is inherently unstable because there is no possibility of stable condition unless the viscosity differences are also very small. These conditions are also shown as pink and brown dotted lines in figure 8.
4. Conclusion
We modelled and mapped the stability dynamics of the VF interaction in dual displacement fronts. Four regimes were identified, and the 
$\tau _{bk}$ can be predicted in each regime. The models and map serve as design tools for harnessing the dual displacement fronts effectively, either for fluid displacement or fluid mixing, depending on its broad applications. We also believe that the models and map can be further developed in the future to explore the instability dynamics of the dual displacement fronts with other configurations. For example, with the different widths of the middle slice, the breakthrough time will be affected due to longer or shorter distances for finger movement and diffusion propagation. Although similar behaviour will still be observed, the regime map will change based on the width slice. Another example is chemical reactions at the interface since they can destabilize the interface or suppress finger development at one of the interfaces. Therefore, this work serves as the first exploration and findings on such an instability dynamics model and map from VF interaction.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.670.
Funding
We acknowledge Japan Society for the Promotion of Science for the funding support with the grant number of 22F22064 and fellowship opportunity for A.P.
Declaration of interests
The authors report no conflict of interest.
