Traveling fronts are recurrent examples of spatio-temporal structures observed in Nature [1–5]. 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 [6–8] 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 [17–19]. 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 [17–19]. In the absence of reaction, spontaneous oscillations can also be observed due to Marangoni instability [20–25]. 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.

The model system consists of a horizontally orientated system (see Figure 1), in which an aqueous solution containing a reactant A, of concentration a ₀, is placed beside an aqueous solution containing a reactant B, of concentration b ₀, along a vertical contact line at time t = 0. The species A and B diffuse at rates D _a and D _b , 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.

The dimensional governing equations for the evolution of the concentrations of the reactants a ̂ , b ̂ and of the product c ̂ 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 v ̂ ̲ = ( u ̂ , w ̂ ) , i.e., ∂ a ̂ ∂ t ̂ + v ̂ ̲ . ∇ ̂ ̲ a ̂ = D a ∇ ̂ 2 a ̂ − k a ̂ b ̂ , (1) ∂ b ̂ ∂ t ̂ + v ̂ ̲ . ∇ ̂ ̲ b ̂ = D b ∇ ̂ 2 b ̂ − k a ̂ b ̂ , (2) ∂ c ̂ ∂ t ̂ + v ̂ ̲ . ∇ ̂ ̲ c ̂ = D c ∇ ̂ 2 c ̂ + k a ̂ b ̂ , (3) ∂ v ̂ ̲ ∂ t ̂ + v ̂ ̲ . ∇ ̂ ̲ v ̂ ̲ = ν ∇ ̂ 2 v ̂ ̲ − 1 ρ 0 ∇ ̂ ̲ p ̂ + g ̲ , (4) d i v ̂ v ̂ ̲ = 0 , (5) or, in dimensionless form, ∂ a ∂ t + v ̲ . ∇ ̲ a = ∇ 2 a − a b , (6) ∂ b ∂ t + v ̲ . ∇ ̲ b = δ b ∇ 2 b − a b , (7) ∂ c ∂ t + v ̲ . ∇ ̲ c = δ c ∇ 2 c + a b , (8) ∂ v ̲ ∂ t + v ̲ . ∇ ̲ v ̲ = S c ∇ 2 v ̲ − ∇ ̲ p , (9) d i v v ̲ = 0 , (10) where δ _b,c = D _b,c /D _a are the two diffusion coefficient ratios, D _c is the diffusion coefficient of species C, and S _c = (ν/D _a ) is the Schmidt number (fixed to 10³ as typical for small species at room temperature in water), with ν = (μ/ρ ₀) the kinematic viscosity, μ the dynamic viscosity and ρ ₀ 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/ka ₀) (with k the reaction rate constant), length L c = D a τ c , velocity U c = ( L c / τ c ) = D a / τ c , concentration a ₀. The pressure is scaled by p _c = (μ/τ _c ) = ρ ₀ S _c D _a /τ _c and a new dimensionless pressure gradient incorporating the hydrostatic pressure gradient is defined as ∇ ̲ p = ∇ ̲ p ̂ / p c − ρ 0 L c g ̲ / p c , with g ̲ = ( 0 , − g ) 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 β = b ₀/a ₀ 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 = ±L _x /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, μ ( ∂ u ̂ / ∂ z ̂ ) = ∂ γ ̂ / ∂ x ̂ at the free surface [21, 27], or in dimensionless form, ∂ u ∂ z = − ∑ i M i ∂ c i ∂ x a t z = L z , (11) where L _x and L _z represent the dimensionless length and height of the system, respectively. In Eq. 11, the dimensionless solutal Marangoni number M _i of species i = (a, b, c) which quantifies the influence of each chemical species on the solution surface tension, is defined as, M i = − 1 μ a 0 D a k ∂ γ ̂ ∂ c ̂ i , (12) where γ ̂ and c ̂ i 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 γ ̂ = γ ̂ 0 + ∑ i ( ∂ γ ̂ / ∂ c ̂ i ) c ̂ i , with γ ̂ 0 the (dimensional) surface tension of the solvent ( c ̂ i = 0 , ∀ i ) . Using Eq. 12, the dimensionless solution surface tension, which is defined as γ = ( γ ̂ − γ ̂ 0 ) / γ ̂ c where γ ̂ c = p c L c , therefore reads γ x , t = − M a a x , L z , t − M b b x , L z , t − M c c x , L z , t . (13)

In Eqs 12, 13, the Marangoni numbers are assumed positive, i.e., M _a,b,c ≥ 0, to describe surfactants decreasing the surface tension of the solvent with ( ∂ γ ̂ / ∂ c ̂ i ) < 0 , ∀ i . From Eq. 13, the surface tension of the initial pure A and B solutions therefore read, γ _A = −M _a and γ _B = −M _b β, respectively.

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

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. M _a = M = M _b , and we assume β > 1 so that γ _B < γ _A . However, the mechanism described below is universal and does not depend on this particular choice.

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 M _c > M _a = M _b , 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 L _z = 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.)

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 (c _S ) and the horizontal component of the velocity field evaluated at the surface (u _S ) (see Figure 4). We note that the profiles of c _S 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 u _S 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 u _S (u _S 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.

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.

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, A _max , 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 A _max in the course of time, the duration of oscillations can also be highlighted (see Figure 6).

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 L _z → 0, we must recover the RD surface tension profiles that cannot oscillate [28]. Hence, there exists a minimum value for the layer thickness L _z above which convection is strong enough to drive the oscillations (numerically calculated to L _z,min = 2 for our specific choice of parameters. Physically, convection is more intense by increasing L _z due to the decrease of the influence of the no-slip boundary condition at the bottom. When increasing L _z above this critical value, the maximum of A _max 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 M _min = 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, L _z > L _z,min , M > M _min , and for intermediate values of δ _c and M _c , where the lower and upper bounds depend on all the model parameters. We have performed the analysis up to (M, M _c ) < 1000, δ _b < 1 and δ _c < 15, L _z < 15, and 1 < β < 15, which defines the range of parameters values tested. We expect that the case β < 1 could be performed similarly.

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 M _a = M = M _b 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. [17–19]. 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.