Exploring buoyancy-driven effects in chemo-hydrodynamic oscillations sustained by bimolecular reactions

Abstract

Exotic dynamics, previously associated only with reactions involving complex kinetics, have been observed even with simple bimolecular reactions A + B → C, when coupled with hydrodynamical flows. Numerical studies in two-dimensional reactors have shown that oscillatory dynamics can emerge from an antagonistic coupling between chemically-driven buoyancy and Marangoni convective flows, induced by changes in density and surface tension, respectively, as the reaction occurs. Here, we investigate reactions increasing both surface tension and density, leading to a cooperative coupling between the flows and show how, in this configuration, buoyancy-driven contribution dampens spatio-temporal oscillations of concentration. We finally identify the key parameters controlling the onset and persistence of the oscillatory instability, namely the density and surface tension gradients, and the systems height.

---

1. Introduction

Chemical oscillations represent captivating examples of a system's spontaneous self-organization and find applications in various fields, including material sciences,^1–6 chemical artificial intelligence,^7–9 environmental science.^10,11 These dynamics are particularly interesting as they also provide understanding of biological functional behaviors^12–24 (e.g. morphogenesis, quorum sensing, calcium signaling, circadian rhythms, synchronization, etc.).

Throughout years, several reactions in which self-sustained dynamics can be induced have been highlighted such as the Belousov–Zhabotinsky,²⁵ Bray–Liebhafsky,^26,27 Briggs–Rauscher,^28,29 Bruk Temkin–Gorodsky Novakovic (BT–GN)^30,31 and chemiluminescence^32–34 reactions. The latter all share a common characteristic essential for inducing oscillatory dynamics: the presence of a nonlinear chemical feedback. Such loops were demonstrated to be crucial in the occurrence of oscillations as in the Brusselator and Oregonator models to name but the most famous.^35,36

However, it was recently shown that even in the absence of complex kinetics (e.g. in A + B → C systems) oscillatory dynamics can occur in spatially extended systems³⁷ when coupled with hydrodynamic effects. The pioneering work of Gàlfi and Ràcz showed that when two pools of reactant solutions A and B initially separated in space react upon diffusive mixing, a reaction–diffusion (RD) front emerges.^38,39 In systems where both reactants present the same initial characteristics (diffusion rate and concentrations), this RD front is localized at the initial position of contact, but as small discrepancies are introduced in the system, the front starts propagating in time.^40,41

In this context, the occurrence of complex dynamics has been investigated in simple bimolecular reactions in the presence of natural convection.^37,42,43 In reactive fluids, local gradients in the physical properties of the system (density and/or surface tension) can be induced as changes in concentration and/or temperature occur in the system, further coupling with the initial RD front and setting the fluid itself in motion. This configuration is modeled in Fig. 1 and will be further discussed in the next section.

Fig. 1 Sketch of the dimensionless A + B → C system. The reactant solutions of A and B, sharing the same properties ([small gamma, Greek, tilde]_A = [small gamma, Greek, tilde]_B = [small gamma, Greek, tilde]_R and [small rho, Greek, tilde]_A = [small rho, Greek, tilde]_B = [small rho, Greek, tilde]_R), are initially separated in space and react upon contact to form the product C in the reactive zone (centered in x₀), increasing locally the surface tension ([small gamma, Greek, tilde]_P > [small gamma, Greek, tilde]_R) and the density ([small rho, Greek, tilde]_P > [small rho, Greek, tilde]_R).

Previous studies have shown that damped oscillations can emerge from a competitive mechanism involving oppositely directed vertical transport processes when the reaction increases surface tension: an upward diffusion relaxation and a downward compressive Marangoni-induced flow. These oscillatory dynamics are further sustained by the addition of so-called antagonistic buoyancy driven flows induced by a local decrease in density and reinforcing the vertical diffusion fluxes [Fig. 2 region I, where ΔM and ΔR are parameters representing the changes in surface tension and density occurring during the reaction, respectively, and will be fully defined in the next section].^37,42,43

Fig. 2 Parametric space diagram (ΔM, ΔR) of the oscillatory behavior in the chemo-Marangoni-buoyancy regime with system dimensions L_z = 20 and L_x = 256 for equal diffusion coefficients of all species and equal properties of all reactants. ΔM and ΔR are parameters representing the changes in surface tension and density occurring during the reaction, respectively. Positive (negative) values of ΔM represent an increase (decrease) of surface tension during the reaction, and positive (negative) values of ΔR represent a decrease (increase) of the density during the reaction. When both parameters are of same sign, the induced flows are acting in opposite directions and the coupling is antagonistic (regions I/III). Contrarily, when ΔM and ΔR are of opposite signs, the coupling is cooperative and the flows reinforce themselves (regions II/IV).

A dynamical regime is defined as oscillatory if the temporal evolution of a variable (concentration, stream function, vorticity,…) presents at least two successive maxima. Oscillations with rapidly decreasing amplitude (of the variable) over time are defined as damped oscillations, whereas oscillations with quasi-constant amplitude over time are defined as sustained. Since the system is closed and will reach equilibrium at very long times, all the oscillations are eventually damped but in the timeframe considered in this study, it is possible to distinguish damped oscillations from sustained ones. If the system is kept out of equilibrium by introducing an inflow of reactants and an outflow of the products, sustained oscillations will persist.³⁷

In the opposite antagonistic case where a decrease in surface tension is coupled to an increase in density [Fig. 2 region III], the competition between both convective flows is essential for the occurrence of oscillatory dynamics.⁴⁴ In both antagonistic cases, sustained oscillations were observed in an optimal range of ΔM and ΔR, eventually dampening and going extinct as one of the convective effects becomes predominant.

Similar effects have been seen, both numerically and experimentally, in non-reactive systems where the introduction of a droplet inducing changes in surface tension led to nonlinear oscillatory dynamics of the surface tension. These oscillations were further affected when buoyancy effects due to a difference in density between the droplet and the solvent were at play, modifying for example the period/amplitude of oscillations or the onset of oscillatory dynamics.^45–47

In this paper, we review the dynamics observed when a local production of C increases both surface tension and density [Fig. 2 region II], completing the classification of oscillatory behavior observed in chemo-Marangoni-buoyancy regime for equal diffusion coefficients and initial properties of all species. Here, both convective flows are acting cooperatively, i.e. we expect the convection rolls to reinforce themselves. In autocatalytic fronts, such a coupling between chemically-induced flows has highlighted particular dynamics (steady asymptotic regimes, also observed in the pure solutal Marangoni effect), different to the ones observed in an antagonistic configuration (spatio-temporal oscillations).^48,49 We note that the opposite cooperative case (i.e. when a local production of C decreases both surface tension and density [Fig. 2 region IV]), has also been investigated but showed no oscillations and will therefore not be discussed.

2. Model

The model system, represented in Fig. 1 is a two-dimensional solution layer of length L_X and height L_Z in the (X, Z) reference frame, where Z is oriented against the gravitational acceleration g = (0, −g). We suppose the reactor isothermal, closed at the bottom and lateral borders and open at the top border with an air–liquid interface supposed non-deformable. Any evaporation is neglected. The system consists of two miscible reactant solutions (A and B) initially separated horizontally as

As the reactant solutions, assumed to have the same properties (density, surface tension and equal initial concentrations), react upon diffusive mixing across the contact line, the product C is formed following the bimolecular A + B → C reaction. The coupling of the latter with molecular diffusion allows for the propagation of a planar reaction–diffusion front, inducing local changes in surface tension and density. The resulting gradients in surface tension and density further trigger convective transport in the system.

The nonlinear dynamics is governed by a set of partial differential reaction–diffusion–convection (RDC) equations for the chemical species A, B, C (eqn (1)–(3)), coupled to the incompressible Navier–Stokes equations for the fluid velocity (eqn (4) and (5))

∂_TA + (V·∇)A = D∇²A − kAB, (1)

∂_TB + (V·∇)B = D∇²B − kAB, (2)

∂_TC + (V·∇)C = D∇²C + kAB, (3)

∇·V = 0, (5)

where T is the time, V = (U, V) is the velocity field, D is the molecular diffusion coefficient assumed constant and equal for all species, P is the dynamic pressure, μ is the constant dynamic viscosity, g is the gravitational acceleration and the symbol 1_Z represents the unity vector in z-direction. The Boussinesq approximation is also introduced⁵⁰ and implies that density changes only affect the gravitational force term . The chemo-hydrodynamical coupling is made via the state equation for the density , and the surface tension, , where I = A, B, C are the dimensional concentrations of the chemical species. In diluted solutions, these equations are assumed to be linear combinations of the chemical concentrations, with and representing respectively the density and surface tension solutal coefficient of the Ith species.⁵¹ρ₀ and γ₀ represent the solvent density and surface tension respectively.

A Marangoni boundary condition (eqn (6)) is imposed at the free surface to describe the chemically-induced shear stress

μ∂_ZU = ∂_Xγ, V = 0 at Z = L_Z. (6)

No-flux boundary conditions are imposed for the chemical concentrations at the four boundaries of the reactor, and no-slip conditions are used for the velocity field at the three solid boundaries.

The system's dimensionless form is expressed by using the reaction–diffusion scales for concentration, A₀, time, t₀ = 1/kA₀, length, , and the derived scales for velocity, , pressure, P₀ = μ/t₀, density, , and surface tension, . The dimensionless solutal Rayleigh, R_i, and Marangoni, M_i, numbers of the Ith species represent the contribution of each species to the density, (∂_Iρ), and surface tension, (∂_Iγ), respectively. They are defined as

using the dimensionless density, , and the dimensionless surface tension, (with i representing the dimensionless concentration of A, B and C).^37,42,43

When the reactants have same initial concentrations and identical diffusion coefficients, the conservation of mass³⁷ implies that a + b + 2c = 1 ∀  x, z, t. This allows to reconstruct the dimensionless concentration field of the product c from the reactants concentrations a and b by introducing two key parameters ΔR and ΔM, defined as, ΔR = R − R_c/2 and ΔM = M − M_c/2 respectively, where R = R_A = R_B, M = M_A = M_B. ΔR (ΔM) tunes the relative importance of buoyancy (surface tension) to the convective flows.

The final dimensionless RDC equations, in which the lowercase characters represent dimensionless variables, therefore read

∂_ta + (v·∇)a = ∇²a − ab, (9)

∂_tb + (v·∇)b = ∇²b − ab, (10)

∂_tv + (v·∇)v = S_c(−∇p + ∇²v − ΔR(∂_xa + ∂_xb)1_z), (11)

∇·v = 0, (12)

with the Marangoni boundary condition

∂_zu = ΔM(∂_xa_s + ∂_xb_s) at z = L_z. (13)

where a_s, b_s represent the dimensionless concentrations at the surface of the reactant solutions. S_c = μ/Dρ₀ = ν/D is the Schmidt number where ν = μ/ρ₀ is the kinematic viscosity. Based on the typical values for the kinematic viscosity ν = 0.0089 cm² s⁻¹ and diffusivity of chemical species D ∼ 10⁻⁵ cm² s⁻¹ in aqueous solutions, S_c is set equal to 1000.

Eqn (9)–(12) are solved numerically by using the alternating direction implicit (ADI) method^52,53 with the defined initial and boundary conditions. Typical simulations are run over a spatial domain of dimensionless length L_x = 256 and variable dimensionless height L_z, discretized over a grid with the integration space steps h_x = h_z = 0.25 and the integration time step h_t = 10⁻⁵.

3. Results

3.1. Explaining the role of buoyancy in the cooperative region II

In the pure chemo-Marangoni regime, when the reaction increases surface tension, the mechanism at the core of the oscillatory dynamics is the competition between the Marangoni-induced flow and the vertical diffusion relaxation. The addition of an upward buoyancy-driven flow (ΔR > 0) reinforces the mechanism leading to sustained oscillations for a range of optimal values of ΔR. In contrast with region I, the addition of a downward buoyancy-driven flow (ΔR < 0) not only fails in providing sustained oscillations in the cooperative case (region II) but also can quench the dynamics. The underlying role of ΔR in this cooperative coupling is described below.

As the reaction occurs, density and surface tension locally increase in the reactive zone, giving rise to a velocity field initially characterized by two convection rolls, one on either side of the reactive zone. The position of the center of a convective roll along the z-axis (z_r) is computed as the position in the system where both the horizontal and vertical velocities are null. In the absence of gravitational forces [Fig. 3a], the convective rolls are centered around 2L_z/3 [Fig. 4 ΔM = 250, ΔR = 0.00] and oscillations develop according to the mechanism summarized above. By contrast, in the other limit case with no Marangoni contribution (i.e. in systems closed at the top surface), the convective rolls resulting from the density gradient are centered at L_z/2 [Fig. 4 ΔM = 0, ΔR = −0.50] but no oscillatory dynamics have been observed.

Fig. 3 Typical spatio-temporal evolution of the product concentration with increasing dimensionless time for systems of length L_x = 256, height L_z = 20 in a (a) pure Marangoni case (ΔM = 300 and ΔR = 0.00) and (b) coupled Marangoni-buoyancy cooperative case (ΔM = 300 and ΔR = −0.50). The fluid velocity field is superimposed on a 2D field of the product concentration c ranging from 0 (blue areas) to the highest concentration, c = 0.44 (red areas).

Fig. 4 Dimensionless temporal evolution of the position of the convective rolls center along the z-axis (z_r) for systems of length L_x = 256, height L_z = 20 and various combinations of ΔM and ΔR. The orange curve represents a pure chemo-Marangoni regime, the green and blue curves represent cases of region II and the red curve represents a pure chemo-buoyancy regime.

When both effects are coupled, as |ΔR| increases the convective cells initially centered at 2L_z/3 are gradually pushed toward the bottom of the system [Fig. 4 ΔM = 250, ΔR = −0.50 and ΔM = 250, ΔR = −1.00]. The amount of product c reaching the surface of the system thus decreases [as can be seen in Fig. 3 at t = 110–150–180 by comparing the convection roll's height in both cases] and, past a certain intensity of the buoyancy-driven flow, becomes insufficient to restore surface tension gradients, further disabling the oscillations. Thus, the coupling with an downward buoyancy flow does not sustain the oscillatory dynamics. Indeed, it dampens them by modifying the velocity field and preventing the renewal of the necessary conditions (i.e. sufficient surface tension gradients) for the emergence of new convective rolls.

The temporal evolution of the local concentration of C is illustrated in Fig. 5 for different combinations of Marangoni and buoyancy contributions. In the pure chemo-Marangoni-driven dynamics (ΔM = 250, ΔR = 0.00), damped oscillations are observed as long as a sufficient gradient in surface tension (ΔM >ΔM_crit) is created by the reaction.^37,42 The addition of antagonistically oriented buoyant flows (ΔM = 250, ΔR = 1.00), allows for sustained oscillations^37,42,43 to be reached in particular ranges of ΔM and ΔR. In the presently studied case, even if both flows reinforce each other, damped oscillatory dynamics are observed for specific combinations of the hydrodynamic parameters. However, the oscillation's amplitude decreases as buoyancy-induced flows are added (ΔM = 250, ΔR = −0.50), further quenching the oscillatory dynamics (ΔM = 250, ΔR = −1.00) due to the pushing of the convective cells to the bottom of the system.

Fig. 5 Dimensionless temporal evolution of the concentration c at a representative point of the system (x₀ − 30, 3L_z/4). The dashed curve represents a typical system in the pure chemo-Marangoni system, the dotted curve represents a typical sustained oscillatory system in region I, the continuous curve represents a typical damped oscillatory system in the cooperative coupling (region II) and the squared curve represent a typical non-oscillatory system in region II.

3.2. Conditions to induce oscillations

The complete picture of the dynamical regimes obtained in region II by varying the convective flows intensity is represented in Fig. 6 and highlights major differences with regions I and III. First, the addition of buoyancy flows allows for the emergence of sustained oscillations in previously studied regions I and III, but can only lead to the dampening of the dynamics in this cooperative case. This shows that the coupling can either have a positive/negative impact on the dynamics by sustaining/quenching them, depending on the changes in density occurring during the reaction. More importantly, it emphasizes that no sustained oscillations are observed in the cooperative case. Secondly, there is a threshold separating non-oscillatory and oscillatory scenarios which scales linearly with the two hydrodynamic parameters ΔM and ΔR.

Fig. 6 Classification of the dynamics observed via a cooperative coupling between chemically induced Marangoni and buoyancy convection in a (ΔM, ΔR) parametric space diagram with L_z = 20 and L_x = 256. The orange squares represent damped oscillatory dynamics and red diamonds represent non-oscillatory systems. The purple line represents the linear threshold between oscillatory and non-oscillatory dynamics.

To observe oscillatory dynamics in the case of a cooperative coupling between flows, several conditions have to be fulfilled. Firstly, the surface tension gradient induced by the reaction must be sufficient to initiate an oscillatory motion, which means that the convection rolls formed must have sufficient horizontal surface velocity to carry a sufficiently large quantity of product away from the center of the reaction front (by the return flow). In particular, we found that at the surface, in both the pure chemo-Marangoni and in the chemo-Marangoni-buoyancy driven systems, the horizontal velocity associated with the first convective roll u^R_s must be greater than the initial horizontal velocity u^I_s (close to the center of the system, in the reactive zone) [Fig. 7]. Since the role of buoyancy is to displace the center of the convective rolls to the lower part of the system, it effectively reduces the horizontal velocity, thereby decreasing the amount of product c brought to the surface and quenching the oscillations.

Fig. 7 (top) Profiles of horizontal velocities at the surface (u_s) at t = 60 for different values of ΔR at constant ΔM = 250. The profiles are displayed in a semi-system (from −L_x/2 to x₀ = 0) since the dynamics is symmetric with respect to x₀ = 0, u(x) = −u(−x). The continuous curves represent oscillatory systems and the dashed curves represent non-oscillatory systems. (bottom) Zoom in the region of interest indicated by the red square. The purple line represent the initial maximum horizontal velocity at the surface (u^I_s) and sketches an approximate limit between the oscillatory and non-oscillatory regimes. u^R_s represents the maximum horizontal velocity associated to the first convective roll. Its value varies over time and is shown here at t = 60.

Based on the Marangoni boundary condition (eqn (13)), it is possible to identify a scaling between ΔM and the horizontal surface velocity. In fact, the typical form of the horizontal component of the Marangoni return flow scales as and fits well with the velocity profiles developed in pure chemo-Marangoni systems.⁴²Fig. 8 shows that the maximum horizontal surface velocity linearly increases with ΔM in the different regimes, as predicted by the analytical solution for the Marangoni return flow u_Ma. A higher slope is observed in the presence of oscillatory dynamics and the velocity field is more intense than in their absence, revealing two distinct linear relations.

Fig. 8 Maximum horizontal velocity at the surface associated to the first convective roll (u^R_s) reached over time (global) as a function of ΔM for different values of ΔR. The velocities are computed in a semi-system (from −L_x/2 to x₀ = 0). The empty squares represent non-oscillatory systems, filled circles represent oscillatory dynamics and dotted lines represent linear regressions.

In this regime II, oscillations are associated with high horizontal velocities as long as there is a dominating effect of the Marangoni contribution (ΔM). High flows do not trivially imply oscillations, as can be seen in Fig. 8 where non-oscillating systems (e.g. ΔM = 450 ΔR = −2.00) present higher maximum horizontal velocities than oscillating systems (e.g. ΔM = 175 ΔR = −0.25). A similar observation applies to previously studied regimes as well.^37,42–44

A second condition arising from the first one is that the product must arrive at the surface in sufficient quantity to recreate gradients and thus re-initiate oscillations. Fig. 9 shows that the maximum vertical velocity v^R in the system decreases as density effects are increased, whatever the values of ΔM considered. The denser the species produced, the lower the maximum vertical velocity, due to the intensification of downward buoyancy-driven flow. As the intensity of the return flow decreases, less product is brought back to the surface of the system, inhibiting the reconstruction of a gradient at the surface and driving the dampening of the oscillatory dynamics.

Fig. 9 Maximum vertical velocity reached over time v^R as a function of ΔM for different values of ΔR in the oscillatory regime. The velocities are computed in a semi-system (from −L_x/2 to x₀ = 0). The dashed and dotted curves represent linear regressions.

Furthermore, Fig. 9 highlights two distinct linear relationships between the flow intensity and the control hydrodynamic parameter ΔM, with the maximum vertical velocity increasing with the surface tension variation. In the pure chemo-Marangoni case (ΔR = 0), a linear relationship is observed between v^R and ΔM. As buoyancy-driven effects are added, a deviation from linearity appears for oscillating systems close to the threshold separating the dynamics (defined by the purple line in Fig. 6) and two linear regimes are thus distinguished. The first one (dotted curves) is close to the boundary where oscillations start to appear and hydrodynamic parameters counterbalance each other. In the second regime (dashed curves), systems are assimilated to pure chemo-Marangoni cases, since density effects are significantly weaker than surface tension effects and are therefore annihilated.

3.3. Influence of L_z

As the amount of product c reaching the surface appears to be a key ingredient for the onset and the control of the oscillatory dynamics, and as it is directly related to the system's height (L_z), we expect L_z to be an important control parameter. We therefore explore the influence of the latter.

Previous studies in the pure chemo-Marangoni case have shown that, in thin layers, the initial convective rolls responsible for the oscillatory dynamics cannot develop and a minimum cell's height (L^min_z) below which wave formation is hindered has been highlighted. Furthermore, L^min_z varies with the initial change in surface tension necessary to induce oscillatory dynamics and scales as L^min_z ∼ ΔM⁻¹. An additional scaling between the period τ at which oscillations are emitted from the center and L_z has been identified⁴² as τ ∼ L_z².

In region I, buoyancy contributions to the flow are at play and the velocity field consists not only of two counter-rotating rolls due to the converging Marangoni flow, but also of two counter-rotating rolls at the bottom of the reactor, pushing the fluid upwardly at x₀.^37,42 Increasing the thickness of the cell, increases the size and intensity of buoyancy-driven forces, which can enhance the oscillatory mechanism (reinforcing the relaxation phase of the cycle) when moderate but, eventually, suppress oscillations as soon as they become dominant and prevail over Marangoni forces. Correspondingly, the characteristic period of the oscillation, τ, first increases as the cell enlarges due to a larger displacement of the concentration field during each spatio-temporal oscillation. However, above a threshold height L^Max_z, linked to a maximum value of τ, the period further decreases linearly as the system gets thicker, marking a switch from Marangoni- to buoyancy-controlled dynamics.

In the case studied here (region II), a range of L_z in which oscillatory dynamics are observed has been highlighted for different values of ΔM and ΔR (Fig. 10). This range increases (decreases) as ΔM (ΔR) increases at constant ΔR (ΔM), between L^Min_z and L^Max_z. In agreement with previous observations, if the system is too thin, the initial Marangoni convective rolls cannot develop, whereas if the system is too thick, the latter rolls are shrunk due to the extend of buoyancy convective rolls and the oscillatory dynamics are inhibited.

Fig. 10 Classification of dynamical regimes in the parametric space diagram (ΔM, L_z) with ΔR = −0.50 and ΔR= −1.50 and L_x = 256. The red diamonds correspond to non-oscillating systems and the orange squares correspond to systems presenting damped oscillations.

However, Fig. 11 shows that the relation τ(L_z) takes the form τ ∼ L_z² for all values of ΔM and ΔR tested, and an increase in ΔM (|ΔR|) translates in an increase (decrease) in the oscillatory period. In contrast to the antagonistic coupling in region I, the period evolves here similarly to the pure Marangoni case. This observation is in line with results obtained in the previous section, showing that the cooperatively coupled cases actually behave like pure Maragoni-driven systems, in which the addition of buoyancy-driven flows has no other effect than to inhibit the complex dynamics.

Fig. 11 Characterization of the oscillation period, τ, as a function of the system height, L_z, for different values of ΔM and ΔR. The solid lines represent linear regressions.

4. Conclusions

In this paper, we have expanded the realm of complex dynamics in bimolecular reactions where Marangoni- and buoyancy-driven flows are triggered by the chemical process. In particular, we have shown that chemo-hydrodynamic oscillations, known to develop in the presence of an antagonistic coupling between Marangoni- and buoyancy-driven contributions, can persist even in a cooperative configuration though the buoyancy effect plays here a quenching role.

The dynamics observed in regions I (antagonistic) and II (cooperative) present morphological similarities, as the core mechanism responsible for the onset of oscillatory dynamics is the same (i.e. introduced in the pure Marangoni case). However, when a bimolecular reaction increases surface tension, the precise nature of the dynamics is determined by the way Marangoni-driven flows are coupled to density effects. When antagonistically coupled (ΔR > 0), the initial oscillatory dynamics can be sustained. When a cooperative coupling is at play (ΔR < 0), the oscillations will always be dampened as the reinforcement of the initial flow pushes the convective cells to the bottom of the system.

The key role of the surface tension gradient in inducing oscillations has already been demonstrated from previous studies.^37,42,44 Here, we deepen the physical understanding of the gradients' key role by showing the importance of the horizontal velocity at the surface associated to the convective roll (u^R_s). In fact, in order for the system to present oscillations, the latter must be greater than the initial horizontal velocity (u^I_s) close to the system's center. Furthermore, distinct linear relations between u^R_s and ΔM have been highlighted in both oscillatory and non-oscillatory regimes. These relations match the analytical scaling predicted for the Marangoni return flow.

In addition, the maximum vertical velocity in the system decreases as the density effects intensify. In turn, less product is brought back to the surface, thus preventing the renewal of a gradient at the surface, which is essential for further oscillations. Moreover, two linear relationships between the maximum vertical velocity and ΔM have been demonstrated, separating two distinct regimes. In the first regime hydrodynamic parameters counterbalance each other, while in the second one coupled systems can be assimilated to pure chemo-Marangoni systems when surface tension effects are significantly stronger than density effects.

Finally, the variation of the system's height L_z shows that oscillatory dynamics are observed only in a range of L_z, which depends on both hydrodynamic parameters ΔM and ΔR. A scaling between the period of oscillation τ and L_z has been identified and takes the form τ ∼ L_z², characteristic of pure Marangoni cases. The latter observation reinforces the idea that the cooperatively coupled cases behaves like pure Marangoni systems.

This work completes the general classification of oscillatory dynamics in a (ΔM, ΔR) plane, paving the way for further theoretical studies in asymmetric systems (different physical properties/diffusion coefficients/concentrations for all the chemical species). This complete picture also extends the range of chemical reactions that can be explored to demonstrate chemo-hydrodynamic oscillations experimentally.

Data availability

The data supporting the findings of this study have been included as part of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Fundings by the PDR-FNRS DYNAMO project and Fondation ULB are gratefully acknowledged. M. A. B. gratefully acknowledges funding from the University of Sassari (Fondo di Ateneo per la Ricerca 2020).

Notes and references

1. G. P. Misra and R. A. Siegel, J. Controlled Release, 2002, 79, 293–297.
2. R. Yoshida, Adv. Mater., 2010, 22, 3463–3483.
3. V. V. Yashin, O. Kuksenok and A. C. Balazs, Prog. Polym. Sci., 2010, 35, 155–173.
4. A. Isakova and K. Novakovic, J. Mater. Chem. B, 2018, 6, 5003–5010.
5. T. Geher-Herczegh, Z. Wang, T. Masuda, R. Yoshida, N. Vasudevan and Y. Hayashi, Macromolecules, 2021, 54, 6430–6439.
6. S. Li, M. M. Lerch, J. T. Waters, B. Deng, R. S. Martens, Y. Yao, D. Y. Kim, K. Bertoldi, A. Grinthal, A. C. Balazs and J. Aizenberg, Nature, 2022, 605, 76–83.
7. D. Przyczyna, P. Zawal, T. Mazur, P. L. Gentili and K. Szacilowski, Jpn. J. Appl. Phys., 2020, 59, 050504.
8. P. L. Gentili, M. S. Giubila and B. M. Heron, ChemPhysChem, 2017, 18, 1831–1841.
9. P. L. Gentili, M. S. Giubila, R. Germani and B. M. Heron, Dyes Pigm., 2018, 156, 149–159.
10. Z. Neufeld, P. H. Haynes, V. Garçon and J. Sudre, Geophys. Res. Lett., 2002, 29, 1534–1538.
11. E. Hernández-Garcia and C. López, Ecol. Complex., 2004, 1, 253–259.
12. A. M. Turing, Phil. Trans. R. Soc. Lond. B, 1952, 237, 37–72.
13. P. Richard, B. M. Bakker, B. Teusink, K. Van Dam and H. V. Westerhoff, Eur. J. Biochem., 1996, 235, 238–241.
14. A. Goldbeter and M. J. Berridge, Biochemical Oscillations and Cellular Rhythms: The Molecular Bases of Periodic and Chaotic Behaviour, Cambridge University Press, 1996.
15. J. D. Murray, Mathematical Biology: I. An Introduction, Springer, New York, 2002, vol. 17.
16. D. Gonze, J. Halloy and A. Goldbeter, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 673–678.
17. S. Jelić, Ž. Cupić and L. Kolar-Anić, Math. Biosci., 2005, 197, 173–187.
18. A. F. Taylor, M. R. Tinsley, F. Wang, Z. Huang and K. Showalter, Science, 2009, 323, 614–617.
19. M. Toiya, H. O. González-Ochoa, V. K. Vanag, S. Fraden and I. R. Epstein, J. Phys. Chem. Lett., 2010, 1, 1241–1246.
20. A. Goldbeter, C. Gérard, D. Gonze, J.-C. Leloup and G. Dupont, FEBS Lett., 2012, 586, 2955–2965.
21. M. Wickramasinghe and I. Z. Kiss, PLoS One, 2013, 8, e80586.
22. G. Dupont, M. Falcke, V. Kirk and J. Sneyd, Models of Calcium Signalling, Springer International Publishing, Cham, 2016, vol. 43.
23. M. A. Budroni, K. Torbensen, S. Ristori, A. Abou-Hassan and F. Rossi, J. Phys. Chem. Lett., 2020, 11, 2014–2020.
24. M. A. Budroni, G. Pagano, D. Conte, B. Paternoster, R. D’ambrosio, S. Ristori, A. Abou-Hassan and F. Rossi, Phys. Chem. Chem. Phys., 2021, 23, 17606–17615.
25. B. Belousov, Sbornik Referatov po Radiatsionni Meditsine, 1958, pp. 145–147.
26. W. Bray and H. Liebhafsky, J. Am. Chem. Soc., 1931, 53, 38–44.
27. L. Negrojević, A. Lončar, J. Maksimović, S. Anić, Ž. Čupić, L. Kolar-Anić and N. Pejić, React. Kinet. Mech. Catal., 2022, 135, 1147–1162.
28. Z. Nagy-Ungvarai, S. Müller and B. Hess, Chem. Phys. Lett., 1989, 156, 433–437.
29. Z. Li, L. Yuan, M. Liu, Z. Cheng, J. Zheng, I. R. Epstein and Q. Gao, J. Chem. Educ., 2021, 98, 665–668.
30. K. Novakovic, L. Bruk and O. Temkin, RSC Adv., 2021, 11, 24336–24344.
31. Ž. Cupić, S. Maćešić, S. Anić, L. Kolar-Anić, A. Ivanović-Šašić and K. Novakovic, React. Kinet. Mech. Catal., 2022, 135, 3–14.
32. K. Pȩkala, R. Jurczakowski, A. Lewera and M. Orlik, J. Phys. Chem. A, 2007, 111, 3439–3442.
33. A. Wiśniewski, M. T. Gorzkowski, K. Pekala and M. Orlik, J. Phys. Chem. A, 2013, 117, 11155–11166.
34. E. Toki, S. Osaki, M. Nakamura, K. Tsuoka, N. Samejima, T. Nakashima, R. Abe, T. Oyama and K. Tsukiyama, Int. J. Chem. Kinet., 2020, 52, 907–917.
35. I. Prigogine and R. Lefever, J. Chem. Phys., 1968, 48, 1695–1700.
36. R. J. Field and R. M. Noyes, J. Chem. Phys., 1974, 60, 1877–1884.
37. M. A. Budroni, V. Upadhyay and L. Rongy, Phys. Rev. Lett., 2019, 122, 244502.
38. L. Gálfi and Z. Rácz, Phys. Rev. A, 1988, 38, 3151–3154.
39. Z. Jiang and C. Ebner, Phys. Rev. A, 1990, 42, 7483–7486.
40. Y. E. L. Koo and R. Kopelman, J. Stat. Phys., 1991, 65, 893–918.
41. H. Taitelbaum, Y.-E. L. Koo, S. Havlin, R. Kopelman and G. H. Weiss, Phys. Rev. A, 1992, 46, 2151–2154.
42. M. A. Budroni, A. Polo, V. Upadhyay, A. Bigaj and L. Rongy, J. Chem. Phys., 2021, 154, 114501.
43. M. A. Budroni, F. Rossi and L. Rongy, ChemSystemsChem, 2021, 4, e202100023.
44. A. Bigaj, M. A. Budroni, D. M. Escala and L. Rongy, Phys. Chem. Chem. Phys., 2023, 25, 11707–11716.
45. N. M. Kovalchuk and D. Vollhardt, J. Phys. Chem. B, 2003, 107, 8439–8447.
46. N. M. Kovalchuk and D. Vollhardt, Phys. Rev. E, 2004, 69, 016307.
47. N. M. Kovalchuk and D. Vollhardt, J. Phys. Chem. B, 2005, 109, 15037–15047.
48. L. Rongy and A. De Wit, J. Chem. Phys., 2006, 124, 164705.
49. M. A. Budroni, L. Rongy and A. De Wit, Phys. Chem. Chem. Phys., 2012, 14, 14619–14629.
50. D. Tritton, Physical Fluid Dynamics, Clarendon Press, 1988.
51. P. Cheng, M. Bestehorn and A. Firoozabadi, Water Resour. Res., 2012, 48, 10902.
52. D. W. Peaceman and H. H. Rachford, J. Soc. Indus. App. Math., 1955, 3, 28–41.
53. C. Fletcher and K. Srinivas, Computational Techniques for Fluid Dynamics 1, Springer, Berlin Heidelberg, 1991.

---