1 Introduction

Traveling fronts are recurrent examples of spatio-temporal structures observed in Nature [15]. In the context of chemistry, fronts are localized reactive interfaces driven by reaction-diffusion (RD) processes that may exhibit unique dynamical behaviors. For instance, for bimolecular fronts, universal time scaling behaviors of the front properties are noted [68] and, for autocatalytic fronts, where a reaction product (the autocatalytic species) catalyzes its own production, spatio-temporal chaos may be observed due to front instabilities [9].

Front propagation is typically influenced by spontaneous hydrodynamic motions, typically buoyancy-driven and/or surface tension-driven (Marangoni) flows, arising from composition (solutal effects) and/or temperature (thermal effects) changes across the front [10, 11]. When such flows are driven across chemical fronts in horizontally oriented-systems (e.g., in Petri dishes or thin layers), many scenarios have been shown to lead to oscillatory phenomena (see Tiani et al. [10] and references therein). In particular, in the presence of Marangoni flows driven across isothermal autocatalytic fronts, transient oscillations of concentration and fluid velocity can be seen before the long-time asymptotic dynamics is reached when the Marangoni number of the autocatalytic species, quantifying its influence on surface tension, is sufficiently large [12]. The antagonistic coupling of surface tension-driven and buoyancy-driven flows can also induce an unsteady periodic behavior of autocatalytic isothermal fronts [13]. Such oscillations, together with oscillations in the temperature field, have also been observed to emerge across exothermal autocatalytic fronts when Marangoni or buoyancy-driven flows are at play due to antagonistic contribution of thermal and solutal effects [13, 14]. The interplay of the oscillating Belousov-Zhabotinky reaction around an A+ B → oscillator configuration with buoyancy forces can also induce new chemo-hydrodynamic instabilities leading to pulsatile fingering and plumes as well as rising or sinking Turing spots [15]. When the same localized oscillating reaction is assumed to change the viscosity of the solutions involved, viscous fingering instabilities can be triggered by the reaction and can induce oscillations in situations that are stable in the absence of chemo-hydrodynamic coupling [16]. While most works on chemical oscillations focused on autocatalytic systems, Budroni et al. recently showed that a transient oscillatory dynamics can be triggered in the presence of Marangoni-driven flows across A+ B → C bimolecular fronts for equal diffusion coefficients and initial concentrations of all species, provided that the surface tension changes during the reaction are large enough [1719]. In that case, such oscillations are explained by the competition of Marangoni convection and the vertical RD relaxation of the front. The combination of both Marangoni and buoyancy-driven flows has further been observed to lead to sustained oscillations around such bimolecular fronts [1719]. In the absence of reaction, spontaneous oscillations can also be observed due to Marangoni instability [2025]. A simple example is the one of a surfactant droplet dissolving under the air/water interface [22]. A solutal Marangoni instability can develop in this system either as steady convection or as oscillations.

In this work, we report an additional mechanism that may lead to transient oscillations and is a unique property of reactive systems where differential diffusion of chemical species combines with chemically-driven Marangoni convection. The mechanism we describe does not require any feedback loop nor prescribed hydrodynamic instability. As detailed below, it is also fundamentally different from the one described by Budroni et al. when the species diffuse at the same rate [17,18], since it does not involve any competition between Marangoni stresses and RD processes as described by the authors and leads to oscillations in the surface tension profiles.

The article is organized as follows. In Section 2, the model system and related governing equations are presented. In Section 3, we describe the emergence of spatio-temporal oscillations of surface tension observed by numerically integrating the governing equations for arbitrary diffusion coefficients and initial concentrations of reactants. Next, in Section 4, we illustrate the control of the oscillatory dynamics as a function of the model parameters. Eventually, conclusions and prospects are drawn in Section 5.

2 Model

The model system consists of a horizontally orientated system (see Figure 1), in which an aqueous solution containing a reactant A, of concentration a0, is placed beside an aqueous solution containing a reactant B, of concentration b0, along a vertical contact line at time t = 0. The species A and B diffuse at rates Da and Db, respectively. The two miscible reacting solutions meet at t > 0 and react according to the isothermal A+ B → C reaction to produce a third species C. The resulting localized region of space where C is produced is called the reaction front. All the species are supposed to affect the surface tension of the solution, therefore inducing gradients of surface tension leading to Marangoni flows. We assume that the air/liquid interface is not deformable and neglect the evaporation processes during the chemical reaction, so that we do not address the dynamics in the air layer. In order to focus on surface tension effects, we also consider the solution density as constant in space and time preventing any buoyancy-driven convection in solution.


FIGURE 1. Sketch of the system. Two aqueous solutions of reactants A and B of concentrations a0 and b0 and of surface tension γA and γB, respectively, are initially separated in space in a 2D domain of length Lx and height Lz. The species diffuse at rates Da and Db.

The dimensional governing equations for the evolution of the concentrations of the reactants


and of the product


in this system are obtained by coupling the reaction-diffusion-convection (RDC) equations for the chemical concentrations to the incompressible Navier-Stokes equations for the dimensional velocity field


, i.e.,


or, in dimensionless form,


where δb,c = Db,c/Da are the two diffusion coefficient ratios, Dc is the diffusion coefficient of species C, and Sc = (ν/Da) is the Schmidt number (fixed to 103 as typical for small species at room temperature in water), with ν = (μ/ρ0) the kinematic viscosity, μ the dynamic viscosity and ρ0 the solution density. To non-dimensionalize the problem, as in Tiani and Rongy [26], we have used the characteristic scales of the reaction-diffusion system: for time, τc = (1/ka0) (with k the reaction rate constant), length


, velocity


, concentration a0. The pressure is scaled by pc = (μ/τc) = ρ0ScDa/τc and a new dimensionless pressure gradient incorporating the hydrostatic pressure gradient is defined as


, with


the gravity acceleration.

The dimensionless initial conditions are separated reactants such that, ∀z, a = 1, b = 0, c = 0, for x < 0 and a = 0, b = β, c = 0, for x ≥ 0, where β = b0/a0 is the ratio between the initial dimensional concentrations of B and A. The dimensionless conditions at boundaries of Figure 1 are no-flux boundary conditions (BCs) for the chemical concentrations at each boundary of the domain. The BCs for the fluid velocity field at the rigid boundaries (x = ±Lx/2 and z = 0) are no-slip conditions, u = 0 = w. At the free surface, we assume w = 0 and we use a Marangoni boundary condition for u derived from the tangential stress balance condition of the form,


at the free surface [21, 27], or in dimensionless form,


where Lx and Lz represent the dimensionless length and height of the system, respectively. In Eq. 11, the dimensionless solutal Marangoni number Mi of species i = (a, b, c) which quantifies the influence of each chemical species on the solution surface tension, is defined as,






are the dimensional solution surface tension and concentration of solute i.

For sufficiently dilute solutions, the solution surface tension is expected to vary linearly with the concentrations. Then, we can write that


, with


the (dimensional) surface tension of the solvent


. Using Eq. 12, the dimensionless solution surface tension, which is defined as




, therefore reads


In Eqs 12, 13, the Marangoni numbers are assumed positive, i.e., Ma,b,c ≥ 0, to describe surfactants decreasing the surface tension of the solvent with


. From Eq. 13, the surface tension of the initial pure A and B solutions therefore read, γA = −Ma and γB = −Mbβ, respectively.

The system dynamics is then obtained by numerically integrating the complete set of Eqs 610 subjected to initial conditions and BCs specified above, with the numerical procedure described in Rongy and De Wit [12]. The length Lx is chosen sufficiently large so that the results are not affected by lateral boundary effects on the time of interest, typically Lx = 350.

3 Spatio-Temporal Oscillations of Surface Tension: Mechanism

When the species diffuse at different rates, oscillations are found in the surface tension profiles for appropriate values of the model parameters (see Figure 2). Such oscillations emerge as local extrema that progressively form at different spatial positions in the surface tension profiles between γA and γB in the course of time. To explain the mechanism leading to such oscillations, we consider the simplest scenario when species A and B have the same influence on the solution surface tension, i.e. Ma = M = Mb, and we assume β > 1 so that γB < γA. However, the mechanism described below is universal and does not depend on this particular choice.


FIGURE 2. Spatial oscillations of surface tension profiles at different times for Ma = 250 = Mb, Mc = 500, β = 2, δb = 0.25, δc = 0.50 and Lz = 10. Such profiles are separated in two plots, (A,B), for clarity.

In the short-time limit, the concentration of C is negligible so that we can neglect its influence on surface tension and the related profiles are nonmonotonic with a global minimum (see Figure 2A at time t = 1). In this limit, from the pure diffusion equations, we can predict the profiles to admit two global extrema (a global maximum and a global minimum) (see reference Trevelyan et al. [28] for the analytical derivation of the corresponding RD profiles, along with the change of notation ρ ↔ − γ). We note that such a global maximum is also formed in the presence of convection but it cannot be seen on the considered scale in Figure 2A due to its negligible amplitude and thus plays no role on the system dynamics. The minimum drives two counterrotating convective rolls, a main convective roll that turns counterclockwise and is of much bigger intensity than the clockwise convective roll on the right (see Figure 3 at time t = 1). Since species A diffuses faster than species B, we note that the main convective roll and surface tension profiles are asymmetric, i.e. more elongated on the side of A than of B. Such an asymmetric convective roll deforms the reaction front across the layer and the maximum production rate of C is found at the surface. As time increases, the effect of species C becomes important. In the reaction zone which is mainly located at the surface, the concentration of C increases. Since Mc > Ma = Mb, the reactants are replaced by the more surface-active product, reducing thereby the surface tension locally and leading to the formation of the first two local extrema (a local minimum and a local maximum) in the surface tension profiles (see Figure 2A at time t = 5). The formation of such extrema locally reduces the intensity of the flow and deforms the inner structure of the main convective roll so that two counterclockwise vortices are formed in the bulk (see Figure 3 at time t = 5). Driven by the vortices, the reactants are transported from their reservoir along the surface at the spatial locations where the flow is the most intense. This generates a new intense reaction zone at the surface (around x = −10 in Figure 3 at time t = 15) while the previous one (located around x = 10 in Figure 3 at time t = 15) reduces in intensity. Inside the newly generated reaction zone, the concentration of C increases, locally reducing the surface tension and thereby inducing two additional local extrema. This RDC process repeats itself progressively forming new intense reaction zones further away towards the left (around x = −20 and x = −30 at times t = 30 and t = 50, respectively) and leading to the observed oscillations in the surface tension profiles. Since the newly generated local minima in such profiles are less pronounced than the previous ones and that the profiles stretch in the course of time, the corresponding vortices are of weaker intensities. Eventually, this process always stops generating new local extrema after some time. In the long-time limit, when all the gradients of surface tension decrease with time, the monotonic properties of the surface tension profiles as predicted by the analysis of RD profiles are recovered Trevelyan et al. [28], here corresponding to two global extrema between γA and γB (e.g., when Lz = 4 and for same values of the other parameters as in Figure 2, the RD profiles properties are recovered numerically after a time of about t = 40. This time typically increases with the intensity of convection as detailed in Section 4.)


FIGURE 3. Focus on the convective rolls centered on the reaction front at different times for the same parameters values as in Figure 2. The fluid velocity field is superimposed on a 2D plot of the production rate which ranges between its maximum value, (ab)max shown in red, and its minimum value, (ab)min = 0, shown in blue. The front is mainly located at the surface and travels in the course of time discontinuously leading to oscillations in the surface tension profiles (cf: Figure 2). The velocity vectors are here tripled compared to their effective length to allow for a better visualization.

Thus, we deduce that oscillations in the surface tension profiles are the results of the subsequent formation of segregated reaction zones along the surface. Each one of them is a source of C production that progressively reduces the surface tension locally as the reaction zones form in the course of time. Such oscillations naturally involve complex spatial dynamics of other observables, such as the product concentration at the surface (cS) and the horizontal component of the velocity field evaluated at the surface (uS) (see Figure 4). We note that the profiles of cS have wave-like tails with multiple maxima located at positions close to the ones of the local minima in the surface tension profiles. For short times, the latter are nonmonotonic with a global minimum and thus, the profiles of uS have a negative and a positive part associated to the counterclockwise and clockwise convective rolls, respectively. As time evolves, the progressive formation of local extrema in the profiles of surface tension come along with oscillations in the profiles of uS (uS indeed oscillates with the gradients of surface tension according to the Marangoni boundary condition, Eq. 12), together with bulk vortices in the system as seen in Figure 3.


FIGURE 4. (A,B) Product concentration cS and (C) horizontal component of the fluid velocity field uS, at the surface, at different times and for the same parameters values as in Figure 2. We note wave-like tails in the profiles of cS and oscillations in the profiles of uS.

We note that the temporal oscillations of surface tension (Figure 5A) and horizontal velocity (Figure 5B) are aperiodic with irregular shapes and amplitudes. Such oscillations dampen in the course of time.


FIGURE 5. (A) Surface tension γ(x, t) and (B) horizontal component of the fluid velocity field u (x, z, t) with z = 0.90Lz, as a function of time at various spatial locations and for the same parameters values as in Figure 2. Aperiodic temporal oscillations emerge and dampen in the course of time with reducing amplitudes as we consider spatial locations further away from the center.

4 Control of the Oscillatory Dynamics

As described in Section 3, the mechanism of oscillating surface tension profiles requires the front to be (mainly) located at the free surface and to be in motion along the surface. As a result, such oscillations are prevented in the symmetric case (i.e., when β = 1 and δb = 1, ∀ M), where two convective rolls of identical intensity leads to a stationary front with surface tension profiles that admit either a global maximum or a global minimum (the symmetric case is described in reference Tiani and Rongy [26]). Hence, oscillations as in Figure 2 can only be triggered in an asymmetric scheme (i.e., when β ≠ 1 and/or δb ≠ 1).

If we arbitrary assume β > 1 (the reverse case β < 1 can be obtained straightforwardly), then γB < γA and the main (biggest) convective roll turns counterclockwise with a surface flow oriented to the left (cf: Figure 3) that brings fresh reactants B towards the side of A. In this case, we numerically find that the condition δb < 1, i.e. that the species A diffuses faster than B, is necessary for the oscillations to occur. Physically, this condition is expected to play the key role of bringing fresh reactants A towards the side of B to feed the reaction zone as it moves along the surface. This explanation is consistent with the variation of oscillations properties (amplitude and duration) with δb as explained below. When β > 1, the condition δb < 1 is therefore found to be essential to initiate and drive the oscillations.

The amplitude and duration of oscillations found in the profiles of surface tension critically depend on the model parameters that modulate both the intensity of convection and the diffusive fluxes of chemical species. To illustrate this, we define the maximum amplitude of the oscillations, Amax, as the maximum difference (in absolute value) of surface tension formed in the positive x-direction between a local minimum followed by a local maximum in the corresponding profiles. By following Amax in the course of time, the duration of oscillations can also be highlighted (see Figure 6).


FIGURE 6. Maximum amplitude Amax of oscillations in the course of time for different values of (A) Lz and (B) δb. Unless varied, the values of the model parameters are M = 250, Mc = 500, δb = 0.25, δc = 0.50, β = 2 and Lz = 6. When increasing Lz or decreasing δb, the maximum amplitude and duration of oscillations in the profiles of surface tension are increased.

The layer thickness is a typical control parameter to modulate the intensity of Marangoni-driven flows and thus, the oscillatory dynamics. When decreasing the layer thickness, convection weakens and so does the driving force for the oscillatory process. In the limit Lz → 0, we must recover the RD surface tension profiles that cannot oscillate [28]. Hence, there exists a minimum value for the layer thickness Lz above which convection is strong enough to drive the oscillations (numerically calculated to Lz,min = 2 for our specific choice of parameters. Physically, convection is more intense by increasing Lz due to the decrease of the influence of the no-slip boundary condition at the bottom. When increasing Lz above this critical value, the maximum of Amax reached in time and the duration of the oscillatory dynamics both increase (see Figure 6A). Similarly, we find a minimum value for M (numerically calculated to Mmin = 200 for the parameters chosen here) since, by increasing M, the difference of surface tension between the solutions of A and B is increased and so is the intensity of convection. Moreover, when varying δb, the relative importance of the diffusive fluxes towards the reaction zone are changed and so are the amplitude and duration of oscillations. In particular, when decreasing δb, more species A diffuse towards the side of B feeding the reaction zone as it moves along the surface. More C is therefore produced and the gradients of surface tension are enlarged. Then, the amplitude and duration of oscillations are found to be larger (see Figure 6B). Oscillations are prevented above a maximum value of δb (numerically evaluated to δb,max = 0.75 for the chosen parameters). Note that the threshold values of the model parameters are provided in this paragraph on an indicative basis only.

Thus, the onset of oscillations critically depend on the model parameters. To summarize, by varying one parameter while keeping the other ones fixed as in Figure 6, if we take β sufficiently large (more precisely, when β > βmin > 1), our numerical results indicate that oscillations emerge when δb < δb,max < 1, Lz > Lz,min, M > Mmin, and for intermediate values of δc and Mc, where the lower and upper bounds depend on all the model parameters. We have performed the analysis up to (M, Mc)


1000, δb < 1 and δc < 15, Lz < 15, and 1 < β < 15, which defines the range of parameters values tested. We expect that the case β < 1 could be performed similarly.

5 Conclusion and Prospects

In this work, we have reported an additional mechanism for the emergence of an oscillatory dynamics unique to reactive systems where differential diffusion effects and chemically-driven Marangoni stresses are at play. Such oscillations are explained by segregated reaction zones that progressively form along the free surface in the course of time. Such segregated reaction zones increase the concentration of C at the surface reducing locally the surface tension leading to the observed oscillations in the profiles of surface tension. This mechanism of oscillations is free of chemical or prescribed hydrodynamic instability.

Next, we have shown the critical influence of the model parameters on the properties (amplitude and duration) and on the onset of oscillations. We have also provided a set of conditions on the model parameters for oscillations to occur with the restriction that Ma = M = Mb assumed for simplicity. A natural extension of this work would then be to relax those assumptions to provide a classification of the oscillatory dynamics in the full parameters space with the motivation to guide future experiments. In this context, while our results could be tested in microgravity or thin film experiments where the effects of buoyancy forces (buoyancy-driven convection) can be neglected [27], we could also include such buoyancy effects into the analysis to extend the application scope of our model to more experimental setups. We could then analyze if the presence of buoyancy forces enhances, reduces or prevents the oscillations of surface tension depending in particular on the Rayleigh number of each chemical species, which quantifies the influence of each of them on the solution density. Also, we have assumed the simplest scenario of an isothermal front traveling in adiabatic conditions (no heat loss to the surrounding). Thermal and/or heat loss effects could also be the subject of future interesting works.

As already mentioned in the introductory part, the spatio-temporal oscillations that arise from differential diffusion effects are fundamentally different from the ones described by Budroni et al. [1719]. In that case, spatio-temporal oscillations in the velocity field and in the concentration profiles are observed when β = 1 and δb = 1 = δc, but there were no spatial oscillation of surface tension (in the sense of the formation and propagation of local extrema in such profiles). Moreover, their mechanism involves an antagonistic effect between the upward relaxation of the front driven by RD processes and the Marangoni downflow that acts to drive the product C in the opposite direction [17,18]. Here, this antagonistic effect is absent. In this sense, our results on the oscillatory dynamics are similar to those found when Marangoni flows are induced across isothermal autocatalytic fronts [12].

We hope that our results will trigger more theoretical and experimental investigations in the growing field of convective effects across traveling fronts.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

LR originally wrote the reaction-diffusion-convection code for autocatalytic fronts [12] which was adapted by RT for bimolecular fronts. The analysis of the numerical results was performed by RT under the supervision of LR. The drafts of this manuscript were written by RT and reviewed by LR. All authors have agreed with the published version of the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer SM declared a past collaboration with the author LR to the handling editor. The handling editor FR declared a past co-authorship with the author LR.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.


The authors thank the “Actions de Recherches Concertées” program and the F.R.S.-FNRS for their financial support.



This article is autogenerated using RSS feeds and has not been created or edited by OA JF.

Click here for Source link (https://www.frontiersin.org/)