A bottom-up approach to construct or deconstruct a fluid instability

Fluid instabilities have been the subject of study for a long time. Despite all the extensive knowledge, they still constitute a serious challenge for many industrial applications. Here, we experimentally consider an interface between two fluids with different viscosities and analyze their relative displacement. We designed the contents of each fluid in such a way that a chemical reaction takes place at the interface and use this reaction to suppress or induce a fingering instability at will. This process describes a road map to control viscous fingering instabilities in more complex systems via interfacial chemical reactions.

A bottom-up approach to construct or deconstruct a fluid instability Darío M. Escala & Alberto P. Muñuzuri * Fluid instabilities have been the subject of study for a long time. Despite all the extensive knowledge, they still constitute a serious challenge for many industrial applications. Here, we experimentally consider an interface between two fluids with different viscosities and analyze their relative displacement. We designed the contents of each fluid in such a way that a chemical reaction takes place at the interface and use this reaction to suppress or induce a fingering instability at will. This process describes a road map to control viscous fingering instabilities in more complex systems via interfacial chemical reactions.
Fluid instabilities at the interface between two fluids are ubiquitous in nature and are responsible for important phenomena that affect everybody's life 1,2 . In some cases, they play a constructive role like in the redistribution of energy in a system but in some other cases, the role is destructive and may pose a serious threat to technical or industrial installations. In most cases, these fluids involve reactants that are known to modify the instability itself [3][4][5][6][7][8][9][10] . Enhanced oil recovery techniques are a clear example where the role of this instability is crucial 11 . A great effort has been done in understanding these phenomena and just recently the effects of reactants have been considered.
We propose in this manuscript a bottom-up approach in order to control instabilities. We focus on instabilities originating at a viscosity jump between two fluids (typically viscous fingering). We develop a system that is likely to produce instabilities and we endow it with the appropriate chemical reactions at the interface that allow us to control the activation or deactivation of the instability at will. In particular, we consider two different fluids with different viscosities and analyze the displacement of one fluid by the other that is being injected into the system. Our results establish the basis to control fluid instabilities that may arise in a broad variety of contexts.
The manuscript is organized as follows. First, we consider the a priori stable configuration with a more viscous fluid displacing a less viscous one and induce an instability at the interface by means of the appropriate chemical reactions. The second part of the paper considers the reverse situation and the interfacial reaction is, then, used to suppress the instability. The final part of the results section is devoted to the full understanding of the mechanism involved and the mathematical model that describes it.

Materials and methods
Experimental section. The system is composed of two solutions, A and B, named displaced and displacing solutions depending on the case of study. One of the solutions, solution A, was prepared by mixing stocks of Poly(Acrylic Acid) (PAA), sodium sulfite, formaldehyde, and pH indicator as follows: The PAA solution was prepared by diluting 1 g of the reagent grade polyacrylic acid with an average molecular weight of 4,000,000 g mol −1 (Sigma), into 180 ml of doubly-distilled water at 80 °C to facilitate solubility. After complete solubilization, the mixture was cooled down to 23 °C and the final volume was kept at 200 ml obtaining a PAA stock solution of 0.5 wt%. The formaldehyde was used directly as a stock solution from a commercial formalin solution (Sigma-Aldrich). The sodium sulfite stock solution was prepared from reagent grade Na 2 SO 3 (Sigma) diluting 25.21 g of the reagent in 100 mL of double distilled water bubbled with argon to avoid oxidation. We included in solution A a pH color indicator (hereafter C.I.) in order to visually monitor changes in pH. We considered as C.I. a 0.4 wt% hydroalcoholic solution of Bromothymol blue prepared by dissolving 1 g of Bromothymol blue sodium salt powder (Sigma) into 50 ml of a 96% ethanol solution diluting up to a final volume of 250 mL by adding 200 mL of doubly-distilled water. The C.I. shows a yellow color for pH values below 6 (acidic state), green color for pH between 6-7 (neutral state), and a blue color for pH values above 7 (basic state) 12 Fig. 1b), the Hele-Shaw cell was made of two square glass plates (25 cm × 25) instead of the plastic plates. This modification was introduced to discard any artifact associated with the cell material and to improve the image quality. To perform the Schlieren measurements, the Hele-Shaw cell was placed between two collimator lenses. The system was illuminated using a pinhole LED light source. An iris cutoff filter was located at the focus point between the CMOS and the collimator lens as indicated in the scheme presented in Fig. 1b 8,18,19 .
Image analysis: circularity. To quantify the effect of the reaction in the obtained patterns, we measure a morphological descriptor called circularity defined as C = 4π[area/perimeter 2 ]. This descriptor indicates how close a region of interest is to a perfect circle. The value of C is calculated directly using FIJI 15 . The methodology is presented in Fig. 2 Fig. 2). Once obtained the binary images, the software calculates the circularity following the definition. As several frames of a complete experimental run are processed, the evolution of the circularity can be plotted as a function of time.
Numerical setup. Simulations are performed using the computational fluid dynamic (CFD) software suite ANSYS Fluent 20.1 20 . The numerical domain is composed of a circular region of radius R = 10 cm discretized using a mapped mesh of radial elements (Fig. 3). The inlet flow velocity at the central boundary is calculated as The injection hole is the circular region of radius r = 2 mm located at the geometric center of the domain. The outlet boundary condition was set as a pressure outlet with p = 0 Pa. The pressure and species fields are discretized using second-order and first-order upwind schemes, respectively. Simulations are performed using the SIMPLE algorithm with a variable time-step setting. The mesh size is chosen after performing a mesh independence study (results are included in the S.I.). The system is simulated using real units to facilitate the comparison with the experiments. The initial conditions are set using the following piecewise functions depending on the experimental case: where N ( r ) is a normally distributed random noise function of amplitude ξ = 0.01 set across the radial coordinate, and ε = 3 mm sets the initial contact region between fluids A and B included for stability reasons.  Table 1.

Results
Direct experiment: viscous solution displaces a less viscous solution. We designed two different solutions with completely different viscosities as described in the Methods section. The large viscosity solution (solution A) contains a large molecular weight (4,000,000 g/mol) Poly(Acrylic Acid) (PAA). The rest of the components of both solutions are selected and distributed in such a way that, at the interface and upon reaction, a dramatic change in the pH takes place. Solution A contains sodium sulfite, formaldehyde, and a pH color indicator as well as PAA, while solution B only contains gluconic acid. In this way, solution A has a significantly larger viscosity 14 : µ A = 150 mPa.s and µ B = 2.2 mPa.s measured at a shear rate of γ r = 50 s −1 . www.nature.com/scientificreports/ In fluid dynamics, the displacement of a fluid by a more viscous one is a stable configuration under any flow rate considered if no reactive processes are involved 3,7,10,26 . This phenomenon is observed in Fig. 4 (lower row labeled Control). Here, solution A (large viscosity) is pumped into the less viscous solution (for this control experiment, solution B is replaced by distilled water to avoid reaction) placed in advance inside a radial Hele-Shaw cell (see Methods for details). Three different snapshots are shown at different moments of the experiment. Note that in all pictures the geometry of the interface is circular, and no instability is observed. The interface becomes blurred as time goes on because both liquids are miscible. Note that in this control experiment, there is no interfacial chemical reaction.   www.nature.com/scientificreports/ Nevertheless, when a change in pH is chemically induced at the interface, the phenomenon is quite different as observed in Fig. 4 (upper row labeled Reactive). Again, the more viscous solution (A) is being injected through the centrally located orifice in the Hele-Shaw cell. Note that almost instantly, the circular geometry of the interface is broken, and some fingers are observed. These fingers grow non-symmetrically until one of them reaches the outer boundary of the cell and the experiment ends. We note a strong dependence of this instability on the inflow rate. As the flow rate of solution A is increased, the instability disappears and a regular circular interface is displayed again. We studied this dependence in Fig. 5 where exactly the same experiment is repeated at different flow rates (only the final configuration of the interface is shown for each case). Larger flow rates produce circular interfaces while for low values of the flow rate the instability appears due to the chemical reaction at the interface.
We conclude that this instability is not a typical fingering instability that it is produced due to the viscous difference at the interface and becomes more unstable as the displacement velocity is increased 10,26 . In our case, the interfacial reaction is crucial to trigger the instability. Apart from the non-reactive control experiment presented in Fig. 4 that shows that the absence of gluconic acid in solution B prevents the instability, several other control experiments were performed and shown in the SI. There, the different reactants were replaced by water and in most cases, part of or the entire interfacial reaction is stopped inhibiting the formation of patterns. Only when the pH color indicator was replaced by water, the instability was observed unaltered (the observation, then, was done using Schlieren techniques described in the Methods section). The conclusion from all these experiments is that the interfacial instability is triggered and maintained by the reactivity at the interface between the two solutions.
A multitude of experiments was performed to characterize this result and they are summarized in Fig. 6. Figure 6a presents the variation of the interface circularity as a function of time for different flow rates. The circularity (calculated as described in the Methods section) provides a value of how close the interface is to a perfect circle (characterized by a circularity value of 1). This value is plotted versus time normalized by the total time of each experiment (time at which the solution A reaches the outer boundary of the Hele-Shaw cell). Note that for low values of the flow rate, the circularity immediately drops to low values. When the flow rate is closer to the transition but still in the instability region, the circularity drops as well but it needs more time to drop. This seems reasonable because as time goes, the average radius of the interface becomes larger and so the linear velocity of the interface becomes smaller. Thus, we move deeper into the instability region. Larger values of the flow rate and clearly out of the instability region are characterized by a constant value of the circularity very close to 1. Figure 6b plots the values of the circularity in the middle of each experiment versus the flow rate. Note that values of the flow rate above Q = 10 mL/min fail to induce instability in the system. Different properties can be measured at the interface to characterize it 4,22,27 . We choose the circularity because it straightforwardly presents the results 8 . A brief description of other properties and some calculations are in the SI.

Mechanism of the induced instability.
We use the Schlieren technique 8,16-18,28 (described in the Methods section) to track changes in the optical index of refraction of the solutions induced by fluid movements that cannot be observed by a direct optical inspection. The use of this technique makes possible the observation of phenomena that are difficult to see with a direct setup. Figure 7 shows pictures of an experiment with a flow rate Q = 3 µL/min. Once the displacing solution is injected (Fig. 7a), both solutions start to react, creating a brownish crust at the interface between both fluids. This crust locally reduces the permeability and increases the pressure inside the cell. This crust produces stagnation areas where more precipitation accumulates. As the displacement solution continues to be injected into de cell, the pressure is increased and, eventually, due to some symmetry breaking mechanism, it is released through some cracks in the crust. This process causes the characteristic patterns observed in the previous figures. The displacing solution contacts, then, fresh displaced solution, and the interfacial reaction starts again, thus, producing more new crust that eventually stops the propagation till the whole process is repeated. This is better appreciated in Fig. 7b in which all the experimental frames between 230 and 730 min were averaged into one single image by image projection methods 15 . This provides dynamic information about the displacing flow in one static image and shows the symmetry-breaking process produced by the accumulation of crust in those regions in which the displacing does not flow. The crust formation is Figure 5. Effect of the flow rate as a trigger for the instability via a reactive interface. All the images are taken at 90% of the final time (t f ) for each case. For Q = 2500 µL/min, the displacing interface remains completely circular. In this case, the characteristic flow time scale is much larger compared to the reaction time scale and a circular stable pattern is obtained. For Q = 200 µL/min, the shape remains mostly circular, but some effects due to reaction and diffusion deformed the initially symmetrical profile. For Q = 5 µL/min, patterns start to be observed forming ramifications around the initially circular shape as described in Fig. 4. For Q = 2.5 µL/min, a fractal-like structure is observed in the displacing solution. Diffusive effects are remarkable for the lower flow rate case. This is observed in the loss of coloration in the displacing solution. The yellowish reaction zone is also observed for all cases, being enlarged for the lower flow rate cases.   www.nature.com/scientificreports/ produced by the agglomeration and subsequent precipitation of PAA molecules driven by the acidic pH of the gluconic acid 29 . In this relatively high acidic condition (pH ≈ 2), the PAA molecules become fully protonated favoring intra-and intermolecular hydrogen bond formation. This phenomenon produces the PAA molecules to aggregate and precipitate almost instantaneously inside the Hele-Shaw cell at the interface. This process was previously reported in Escala et al 8,14 in a system with a similar composition. More information regarding the precipitation formation is included in SI. Figure 7c shows another mechanism that is observed to play a major role. Once the injection is stopped a chemical reaction wave is observed to propagate in the direction opposite to the previous flow velocity. In fact, this chemical wavefront is always present in these experiments and it is responsible for the dynamics of the fingers.
The mechanism of the chemical front needs two conditions to occur. One condition is that SO 3 2− must be in solution A (see Supp. Info. for more details). The acid pH of the gluconic acid solution interacts with the SO 3 2− of solution A and producing HSO 3 − (which is essentially acid) by equilibrium displacement as expressed by the following equations 8,14 : where PA − is the polycarboxylate ion.
It is also important to note that there is a large gradient of protons between solutions A (pH ≈ 12) and B (pH ≈ 2). This explains why the observed chemical front is indeed an acid front, and the color indicator switches from blue (basic state) to yellow (acid state) (see Figs. 5 and S16 in Supp. Info). The second condition for the reaction front to occur is the existence of a large difference in diffusivity between the polymer molecules (D PAA = 1 × 10 -11 m 2 /s) 4,8,21,23 and the H + (D H + = 9.3 × 10 -9 m 2 /s) 24,25 of the solution B. This difference is enhanced by the effect in the proton apparent diffusivity produced by the PAA molecules which act as a reversible proton acceptor 30,31 . This difference in diffusivity plays a major role as a condition for the front propagation when autocatalysis is not possible. As no oxidant species are in the medium, the nature of this chemical front is far from the classical autocatalytic approach vastly observed in many pH-oscillators 30,31 . Diffusion also plays an important role in the system stability 32 (See SI for more details).
From these observations, it seems straightforward that there must be a coupling between the flow velocity and the chemical rate. In order to stress this, we estimated the non-dimensional Damköhler number following 8,22,27 (see Supp. Info for more details). Its physical meaning is just the ratio of the hydrodynamic characteristic time over the temporal scale for the reactive processes involved. Figure 8a plots the estimation of the Damhköler number (Da) as a function of the average radius of the interface for three characteristic flow rates. The average interface radius (R) (Fig. 8c) is defined as the minimum radius at the blue-yellow boundary in the first experiment www.nature.com/scientificreports/ in which pattern formation was observed (Q = 5 μL/min). Note that as the experiment evolves, the interface average radius becomes larger and so the hydrodynamic temporal scale becomes smaller and Da larger. Figure 8b plots all the Damhköler numbers estimated for all the experiments performed calculated for an average interface radius of 20 mm. For Q > 10 µL/min, Da < 1 indicating the prevalence of the advective process, while for Q < 10 µL/min, Da > 1 demonstrating the prevalence of the reactive process. Da > > 1 for Q = 0.5 µL/min and Da < < 1 for Q = 2500 µL/min are the extreme cases that confirm the tendencies. We conclude that for larger injection flow rates, the hydrodynamics is such that the reaction front has no practical effect and the normal behavior is observed (no instability). Only when the reaction is allowed to play a role (Da > 1) the instability is seen. These results demonstrate the physical aspects of pattern generation. They also show that there are two chemical processes involved playing a major role in pattern formation and directly affecting the physics of the system. Once solution A is injected, there is a reaction front that moves opposite to the flow direction. As this is a radial system, the flow velocity decreases with the radius. Once the injection velocity matches the velocity of the chemical front, the crust starts to accumulate and, thus, reduces locally the permeability of the system. This permeability loss locally increases the pressure making the displacing fluid break the crust and be ejected through the cracks generating the characteristic fingers. When both fluids make contact again, the reaction front forces the finger to recoil. The whole process is better appreciated in the movies included as supplemental info.

Reverse experiment: less viscous solution displaces a viscous solution.
We consider in this section the reverse situation. The more viscous solution (solution A) is located initially inside the Hele-Shaw cell and the less viscous solution (solution B) is injected into the cell with a constant flow rate becoming, thus, the displacing solution. This is an unstable situation from a fluidic point of view and in the absence of reaction produces a viscous fingering instability that has been widely characterized in the literature and appears in many industrial applications with negative impact 7,33,34 . We analyze here the role played by the interfacial reaction in controlling the instability. Figure 9 shows the results of such an experiment for two different values of the injection velocity. Figure 9a (lower row) presents the non-reactive control experiment where the injected solution is replaced by distilled water. This is an unstable configuration and viscous fingering is immediately observed. The upper row presents the case where the interfacial reaction is allowed (the displacing solution is gluconic acid as described in the Methods section). In this case, the evolution is completely different, the instability is suppressed and there is an effective displacement of the more viscous solution A. In Fig. 9b, the same experiments are repeated but at larger flow rates. Note that in this case, the instability is not suppressed by the interfacial reaction and it is observed both in the control experiment and in the reactive one. Similar to Fig. 5, the effect of the flow rate is observed in Fig. 10. Here, snapshots taken for different experiments with different flow rates are compared. Experiments were run (as explained in the Methods section) till the displacing solution touches the external wall (labeled t f in the manuscript and figures). In this figure, we compare pictures taken at 0.25 t f . Note that for Q > 10 mL/min the interface becomes unstable. This corresponds with Damköhler numbers below 1 as described in the previous section. We characterize the intensity of the instability by measuring the circularity of the interface and the results are plotted in Fig. 11. Figure 11a presents the variation of the circularity along with the experiment for three different flow rates. For Q = 10 mL/min (red line in Fig. 11a) the value of circularity is high and almost constant during the entire experiment. In the two other cases, the interfacial reactions are not able to suppress the instability, and circularity drops to very low values (close to zero). As time evolves, circularity increases due to a combination of diffusive effects and the chemical front previously described. Note that in this case, the reaction front moves in the same direction as the flow rate. Figure 11b plots the values of the circularity for several flow rates measured at t = 0.5 t f . There is a clear Figure 9. Stabilization of an initially unstable front for (a) Q = 10 mL/min and (b) Q = 200 mL/min. In both flow cases, a reactive (upper row) and a non-reactive (lower row) experiment were conducted. Non-reactive cases (control) are always unstable leading to the development of a fractal pattern that rapidly reaches the border of the observation region. The reactive cases are sensitive to the flow rate Q, in (a) the displacing front growths stable showing a more circular pattern. When the flow rate is increased (b), a fractal shape is observed again. Note that in this case, once the border of the reactor is reached, the fractal ramifications increase their thickness changing the front shape. This is not observed in the control case, where once the pattern reaches the border of the reactor, it remains fractal. www.nature.com/scientificreports/ transition as the flow rate is increased (above 10 mL/min), below this critical flow rate the fingering instability is suppressed while above the instability is still present. Note that this transition is marked by the Damköhler number. Da > 1 corresponds to Q < 10 mL/min and a stable interface is observed while for Da < 1 (corresponding to Q > 10 mL/min) the interface is unstable. Another important feature, especially relevant for industrial applications, is the total amount of displaced solution per experiment. This is shown in Fig. 12 for the control case without reaction and the reactive situation (Q = 10 mL/min in both cases). Note that once the control case, with fingering instability, reaches the cell boundary, solution A stops being displaced (black curve in the figure). On the other hand, when the reactive situation is considered and the instability is suppressed, the displacement of solution A continues for larger times, and the total volume of solution A displaced is significantly larger (an additional 123% more in our case. Red curve in the figure).
Similar to the previous case, we also use the Schlieren view to understand the stabilization mechanism. The stabilization process of an initially unstable configuration is shown in Fig. 13 through the Schlieren technique. Note that the mechanism proposed for the direct experiment above is still valid in the present case. For low flow rates, where chemical and advective timescales are comparable (Da≈1), the chemical front moves synergistically with the flow front. The reaction produces polymer precipitation and the creation of an effective crust. As pressure increases due to the incoming flow, the crust breaks. At this point, the faster reaction rebuilds the crust. This process, repeated randomly around the interface, results in a stable rounded interface propagating smoothly. When the flow rate is increased (and the Damköhler number is smaller than 1), the fluid moves faster than the reaction characteristic time, and no effective wall is created, thus, the system remains unstable (see S.I. for some pictures). Reacting interfaces for several flow rates (Q). All images are taken at t/t f = 0.25 for each case. As the flow rate is decreased, the initially unstable front changes from a fractal-like structure into a more stable circular shape. For intermediate flows (Q = 15 and 50 µL/min) the fractal geometry is still there but the fractal fingers get enlarged due to the reactive process. Only for the lower flow rate case (Q = 10 µL/min), the interface becomes rounded and stable. Figure 11. Quantitative measurements of the morphological changes of the interface as a function of the flow rate Q. (a) Evolution of the circularity for three different flow rates: 10, 100, and 200 µL/min. For Q = 100 and 200 µL/min, the circularity rapidly drops to values near zero for t/t f < 0.25. This is due to the fractality of the displacing front. For t/t f > 0.25, the circularity increases with time due to reactive effects as shown in Fig. 9b. For Q = 10 µL/min, the circularity remains stable between 0.5 and 0.6 during the recorded time. Only some selected cases are plotted here for clarity (see the rest, together with the corresponding control experiments, at the SI). (b) Semi-log plot of the circularity for all the studied cases at t/t f = 0.5. The system becomes more unstable as the flow rate increases. For Q = 10 and 5 µL/min, the system becomes more stable due to the interfacial reaction.

Mathematical model and numerical solutions.
Here, we incorporate the mechanism described above into a mathematical model that couples the dynamics of the displacement of reactive fluids in a porous media inside a Hele-Shaw cell, the existence of a chemical front, and polymer precipitation that locally affects the permeability of the system. We considered parameter values as measured experimentally.
For a system with viscosity µ, porosity ϕ and permeability κ, the governing equations, based on Darcy's Law for incompressible fluids, are given by the following Reaction-Diffusion-Advection set of equations 10 ,  . Front stabilization mechanism observed using the Schlieren technique. Starting from an initially unstable situation, the stability of the system is increased by the synergy between the polymer aggregation (reactive front) and the flow displacement. This is only observable when the reactive and flow fronts move synergistically.
Scientific Reports | (2021) 11:24368 | https://doi.org/10.1038/s41598-021-03676-z www.nature.com/scientificreports/ where � r = x, y is the position vector. We considered a simplified version of the reactive solutions where A stands for solution A containing the polymeric solution, B for the solution B (containing the gluconic acid), and C for the polymer precipitate. Note that species A, B, and C depend both on space location and time, A = A(� r, t) , B = B(� r, t) and C = C(� r, t) . ϕ is the porosity of the system and it is set constant for a Hele-Shaw cell 8,10,26,32 . D A , D B , and D C are the diffusion coefficients of A, B, and C respectively. We assume that the polymeric solution (A) diffuses slower in comparison with the acidic solution (B) and that the precipitate (C) does not diffuse at all 8,23,35,36 .
The permeability of the system κ(C) depends on the polymer precipitate and it is modeled by the following equation 35 , where κ 0 is the permeability when C = 0 (no precipitate in the medium), and it is calculated as κ 0 = h 2 /12 26,[37][38][39] , and h is the separation gap of the Hele-Shaw cell. We define κ m = κ(C m ) as the permeability when C = C m and R κ = ln(M 0 /M m ) as the permeability log-mobility ratio, where M 0 = κ 0 /μ and M m = κ m /μ are the mobilities when C = 0 and C = C m , respectively. C m is a reference value that we choose as the maximum possible concentration of the species C in the medium. The parameter R κ quantifies the influence of precipitation on permeability changes. A positive value of R κ indicates that C reduces the permeability of the porous matrix locally 35 . This is the situation assumed in our case, as the precipitate formation clogs the cell reducing the permeability. For this work, a value of R k = 80 is chosen ad-hoc based on experimental observations. We consider that the viscosity of the system µ(A, B) is only affected by the displacing and displaced solutions and it is defined by the following relation 10 : where μ A = μ(A 0 ,0,0), represents the viscosity of the polymeric solution for a specific concentration A 0 in the absence of any other solution in the medium. The parameter R B compares the viscosity of the two reactant solutions 10 . Following the experiments, both parameters µ A and R B are fixed, and the stability of the system is only affected by the initial configuration. Those values are also estimated from experimental observations.
A simplified model is used to simulate the reactivity between the displacing and displaced solutions. This model is based on the mechanisms proposed in the previous section and it is given by the following set of chemical equations: where Eq. (7) models the reactive front. It represents the bisulfite generation due to the acidity of solution B in one single step. D represents any secondary product species generated by this mechanism and it is included to preserve mass conservation. By sake of simplicity, the concentration of D is set constant and equal to zero during the entire simulation. Equation (8) models the polymer precipitation (C) that is produced when both solutions (A and B) make contact at the interface. As previously mentioned, the polymer precipitation is produced by aggregation of PAA molecules due to a low pH interfacial condition 29 . The use of such an equation to model precipitation has been extensively used in many previous works 10,28,35,40,41 . k 1 and k 2 are the rate constants. k 1 was calculated from the experimental front velocity (see Supp. Info for more details). The value of k 2 was set as k 2 = k 1 /10 to better fit the experimental observations. This set of equations is integrated as described in the Methods Section. Figure 14 presents the numerical results obtained for the direct experiment. Figure 14a shows the overlaid image of the concentration fields of species A and C for Q = 5 µL/min. Note that both the pattern structure and the polymer precipitation are qualitatively well reproduced by the model (in the S.I. some animations of this simulation are shown with a good agreement with the experimental ones). In the non-reactive case (lower row in Fig. 14a, labeled as control), the system remained stable as in the experimental counterpart. Figure 14b shows the circularity evolution in normalized time (t/t f ) for three representative cases (the remaining cases are presented as supplementary information) in good agreement with the experimental results in Fig. 6. When the flow rate is increased, the system reactivity results diminished compared to the advection and the system remains stable (circularity ≈ 1). For lower flow rates, the chemical timescales match the advective characteristic times at some specific radius. This, in combination with the local reduction of permeability (and the subsequent pressure increment) produced by the precipitation, leads to the fingering formation. Note that in the simulations the initial circularity is always close to 1.
The circularity variation as a function of the flow rate is presented in Fig. 14c, also in good agreement with experiments. Figure 15 shows the numerical results for the reverse case. Figure 15a shows a comparison between two different flow rates Q = 10 and 200 µL/min for both reactive (upper row) and non-reactive (lower row) cases in normalized time. Figure 15b shows the circularity variation versus the normalized time for three representative flow rates (the remaining cases are included as SI). Figure 15c shows the circularity variation as a function of the  Figure 16 shows the total volume of the displacing solution that is recovered at the end of each simulation and the values are compared with the control non-reactive experiment (compare with the experimental Fig. 12).

Discussion and conclusions
We considered in this contribution a highly viscous polymeric solution and studied the problem of its displacement by another miscible less-viscous solution. We designed the second solution in such a way that a reaction takes place at the interface that locally changes the stability mechanism. In this way, the stable configuration (viscous solution pushing a less viscous solution) becomes unstable and vice versa.
The physical mechanism underlying is a dramatic decrease in the pH at the interface that is induced once the two solutions interact and only at the interface. This change in pH results in two different chemical processes that combined with the fluid displacement are responsible for the destabilization/stabilization of the system: precipitation and the reactive front.
The precipitation is produced by the alteration of the spatial configuration of the polymer at low pH. This process is driven by the protonation of the PAA molecule that facilitates the aggregation and precipitation by hydrogen bond formation. In the direct experiment, the precipitate accumulates in the region between the two solutions and locally reduces the permeability behaving like a stopper. This, consequently, increases the pressure inside the cell. At a certain level of pressure, the crust breaks and the displacing solution is ejected. This break happens by spontaneous symmetry break. In the reverse case, (the less viscous solution displacing a more viscous one), the whole mechanism is the same, just the polymer aggregation process and the flow displacement act synergistically to produce a stable displacement.
The reaction front is produced by a combination of several reaction-diffusion processes. In the first place, the acidic character of the front is explained by the generation of bisulfite at the interface via equilibrium displacement. This reaction is potentiated by the large gradient of protons between solutions A and B.
The Damköhler number (the ratio of the hydrodynamic characteristic time over the temporal scale for the reactive processes involved) is used as an indicator that marks the region of flow rates that can change the stability of the process. Large Damköhler numbers are dominated by significant reactive processes that modify the stability at the interface as the reaction has time to compete with the hydrodynamic process. We also characterized the role of the diffusion against advection and reaction by studying the Péclet and the Péclet-Damköhler Color codes are as follows: blue shows maximum of A concentration while light green is for the whole domain occupied by B specie, the concentration of the precipitate (overlayered in the same figure) goes from white (no precipitate) to red (maximum of precipitate). All results are presented in normalized time. (b) Circularity variation of the same three representative cases presented in the experimental Fig. 6a as a function of the normalized time (t/t f ). The remaining cases are shown in the SI. (c) Circularity variation of all simulated cases as a function of the flow rate Q at t/t f = 0.5. Simulation parameters are indicated in Table 1 www.nature.com/scientificreports/ numbers, respectively. We demonstrated that differential diffusion is fundamental for the front to occur, and it is mainly produced by the PAA molecule that acts as a proton acceptor. Due to the similarities of the patterns here presented with those observed in dissolution fronts instabilities 42 , other specific numbers, such as Zhao number, may be applied to characterize this system as well [43][44][45] . However, more research is needed in order to adapt the characteristics of an ideal Hele-Shaw cell system into a more specific application. This will be addressed in future work. The mechanism observed in the experiments has been translated into a mathematical model and numerically integrated. The results show an important agreement with the experimental evidence demonstrating the validity of the proposed mechanism.
The mechanism used along this manuscript took advantage of the pH-induced precipitation but eventually and depending on the application, different not-involving-precipitation approaches could be imagined. The protocol proposed in this paper opens a path to design reactions that change the stability properties at the interface by synergistic mechanisms between chemical and advective processes. The possible applications are countless in industrial configurations as well as in processes of recovery of natural resources.  Table 1 Figure 16. Quantitative numerical comparison between the displaced volumes for the reactive and control cases when Q = 10 µL/min. The total displaced volume is compared at the end of the simulation (when all displacing is injected). Note that the total displacement is improved up to 174% comparing with the nonreactive case. The center vertical line indicates the time when the pattern in the non-reactive case reaches the boundary of the domain. The simulation parameters are indicated in Table 1.