# Salinity-driven convection mechanisms in stratified brines

Christian O Oko, Berin Šeta, Aliaksandr Mialdun, Mounir M Bou-Ali, Valentina Shevtsova, Arnault Lassin, Henri Bataller

► To cite this version:

Christian O Oko, Berin Šeta, Aliaksandr Mialdun, Mounir M Bou-Ali, Valentina Shevtsova, et al. Salinity-driven convection mechanisms in stratified brines. Chemical Engineering Science, 2026, 324, pp.123224. 10.1016/j.ces.2025.123224. hal-05460396

HAL Id: hal-05460396

https://brgm.hal.science/hal-05460396v1

Submitted on 15 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Salinity-driven convection mechanisms in stratified brines

Christian O. Oko ID a, Berin Šeta ID b,*, Aliaksandr Mialdun ID c, Mounir M. Bou-Ali ID c,
Valentina Shevtsova ID c, Arnault Lassin ID d, Henri Bataller ID a

## ARTICLE INFO

### Keywords:

Brine

Electrolytes

Diffusion

Interferometry

Double-diffusion convection

Instability

## ABSTRACT

Understanding transport and mixing in stratified saline systems is critical for predicting the behavior of brines in natural aquifers, industrial reservoirs, and engineered disposal sites. These multi-ion solutions often exhibit complex instabilities driven by differential diffusion and compositional gradients. The onset and morphology of such convective mixing remain poorly predicted. We investigated double-diffusive (DD) and diffusion-layer convection (DLC) in superimposed aqueous solutions of the salts typical of saline aquifers, sodium chloride (NaCl) and sodium sulphate (Na₂SO₄). The study combines thermodynamic modeling, optical interferometry experiments, and nonlinear numerical simulations to explore convective instabilities in a ternary system. Our findings reveal a rich variety of convective scenarios depending on salt configuration and concentration ratios. When the faster-diffusing NaCl was placed above Na₂SO₄, diffusion-layer convection occurred with a delayed and asymmetric onset of instability, an experimentally demonstrated feature not reported previously. In contrast, when the positions were reversed, the system developed double-diffusive fingers that grew slowly due to the small density ratio. These fingers exhibited an unusual morphology, consisting of extremely fine, vertically textured structures that gradually merged away from the interface. This formed a large area of diffuse mixing and suppression of coherent convective structures. In all cases, classical stability criteria failed to fully predict the onset and nature of convection. Instead, we identified the initial position of the system on the stability map, as determined by the full diffusion matrix, as a critical factor.

## 1. Introduction

Diffusion controls mass transport in a wide variety of geochemical, environmental, and industrial processes involving aqueous electrolytes. These include sediment diagenesis, saline lake evolution, mineral growth, and, of increasing relevance, the geological storage of CO₂ in saline aquifers, where induced salt precipitation hinges on a precise understanding of diffusive transport in concentrated brines (Applin and Lasaga, 1984; Steefel and Lichtner, 1994; Yadav and Kim, 2015; Luo et al., 2022; Yan et al., 2025; Rehman et al., 2022).

Sodium chloride (NaCl) and sodium sulfate (Na₂SO₄) are major constituents of many natural brines, including those targeted for geological CO₂ storage. Consequently, their thermodynamic properties and binary diffusion coefficients in aqueous solutions are well-characterized (Archer, 1992; Pitzer, 1991; Hall et al., 1953; Rard and Miller, 1979; Rard, 2021). However, predicting diffusive transport in their ternary mixture, NaCl + Na₂SO₄ + H₂O, is critical for assessing brine evolution

and salt precipitation during CO₂ injection, yet it remains a significant challenge. Although diffusion coefficients have been reported for numerous compositions (Rard et al., 1996; Albright et al., 1998; Annunziata et al., 2000; Fu et al., 2001, 2002), critical gaps persist: (1) data are lacking for the specific compositions relevant to this study, and (2) established predictive models are inadequate. The Nernst-Hartley formalism is the standard model for predicting ternary diffusion coefficients. It calculates these coefficients from ionic conductivities at infinite dilution, applying a viscosity correction for high concentrations. However, this approach has been shown to fail for the NaCl + Na₂SO₄ + H₂O system (Annunziata et al., 2000; Fu et al., 2001).

To address these gaps, we developed an integrated approach combining a new thermodynamic model with direct experimental validation. First, we implemented a PHREEQC-based model to calculate diffusion coefficients for the NaCl + Na₂SO₄ + H₂O system, building on our previous work for binary solutions (Oko et al., 2024). This model provides an alternative to the inadequate Nernst-Hartley formalism for the

<footnote>a Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LFCR UMR5150, Anglet, France</footnote>

<footnote>b Dept. Civil and Mechanical Engineering, Technical University of Denmark, 2800, Lyngby, Denmark</footnote>

<footnote>c Fluid Mechanics Group, Faculty of Engineering, Mondragon University, 20500, Mondragon, Spain</footnote>

<footnote>d BRGM, F-45060, Orleans, France</footnote>

<footnote>* Corresponding author.
E-mail address: berse@dtu.dk (B. Šeta).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic illustration of the experimental setup. a) The diffusion cell with the central square indicating the field of view. b) A Mach-Zehnder interferometer. This setup enables the measurement of interference fringe patterns (interferograms) caused by the superposition of a reference beam and a test beam which passes through the cell.</figcaption></figure>

compositions of interest. Second, to validate the model, we performed double-diffusion experiments in a Hele-Shaw cell, using optical interferometry to observe the mixing between gravitationally stable layers of NaCl and Na₂SO₄ solutions. The instabilities observed in these experiments serve as a sensitive benchmark for our numerical simulations, which were conducted using the calculated diffusion coefficients.

By combining a robust thermodynamic model, a numerical simulation framework, and direct experimental validation for targeted compositions, this work extends the understanding of diffusion and instabilities in the complex, environmentally relevant NaCl + Na₂SO₄ + H₂O system.

## 2. Materials and methods

In this work, we examine two characteristic salts of saline aquifers. A layer of sodium chloride (NaCl) aqueous solution is superimposed on a layer of sodium sulfate (Na₂SO₄) aqueous solution or vice versa, varying both the concentrations and their ratios.

### 2.1. Solutions preparation

For all the experiments in this study, NaCl (Sigma-Aldrich, 746398-500G, anhydrous, ACS reagent, purity 99%) and Na₂SO₄ (Sigma-Aldrich, 239313-500G, anhydrous, ACS reagent, purity 99%) were used without further purification. De-ionized water retrieved from a Millipore Milli-Q filtration station was degassed and used to prepare a volume of 100 mL for every saline solution to the required concentration in mol/L. The required salt masses to prepare the solutions were weighed using an analytical balance (Sartorius TE313S, with a resolution of 10⁻²g/200 g). All solutions were prepared and experiments were performed at atmospheric pressure and a room temperature of T = (22.0 ± 0.5)°C, which was maintained with the help of a dedicated air conditioning system in the laboratory.

### 2.2. Experimental methodology

We employed interferometry, a non-invasive optical technique, to monitor subtle concentration changes by detecting variations in the refractive index.

A sample cell was placed in a Mach-Zehnder interferometer (Fig. 1), where interference patterns were generated by splitting and recombining a coherent light beam. This allowed for precise mapping of the refractive index field over time, providing high-resolution insights into concentration dynamics. A polarized He-Ne laser beam (λ = 632.8 nm, 5.0 mW) is enhanced, collimated, and expanded by a spatial filter. The collimated beam is further expanded by two lenses, and then split into two beams by a beam splitter. The reference beam passes through free

ambient air, while the test beam (measurement beam) passes through the diffusion cell. The beams are recombined at a second beam splitter to form an interference pattern, which is generated by introducing a phase shift achieved by tilting mirrors 1 or 2. The pattern is magnified by a convex lens and captured by a CCD PIKE camera (Allied Vision Technologies, PIKE-F421B), capable of recording at 15 frames per second with a full frame resolution of 2048 × 2048 pixels.

The diffusion cell is a double-channel flow-through parallelepiped quartz cuvette with internal length Lₓ = 17 mm and height Lᵧ = 58 mm, with an optical path l = 1.05 mm. The custom-designed cell was ordered from Aireka Scientific Co. Ltd. (iCuvets.com). It was filled in two steps. First, the cell was completely filled with a less dense fluid through the bottom inlet. Second, the cell was simultaneously filled with the less dense solution from the top inlet and the denser solution from the bottom inlet, displacing the initial solution and ensuring that the low-density fluid remained on top of the high-density one. Once filling was complete, we started recording images.

Image acquisition is managed by LabVIEW software developed in-house. Typical interferograms are shown in Fig. 2(a1) and (a2). In our experiments, the acquisition frequency is 1 fps (1 Hz). The field of view (the central square indicated in Fig. 1(a) includes 2048 pxl × 2048 pxl or 10.936 mm × 10.936 mm. Thus, the spatial resolution is (5.34 ± 0.05)µm per pixel. The acquisition rate was intentionally set to 1 fps to match the characteristic timescale of the physical processes under investigation, which evolve on the order of several seconds to minutes. Although the camera is capable of higher frame rates, increasing the acquisition frequency does not provide additional relevant physical information for the present experiments. Moreover, a rate of 1 fps ensures stable real-time image processing without generating unnecessary data redundancy or overloading the integrity of the automated protocol. This choice represents a compromise between temporal resolution and data handling.

The recorded interferograms are post-processed (see Fig. 2) by first performing a classical 2D Fast Fourier Transform (2D-FFT) using MATLAB algorithms to convert the images from the spatial domain to the Fourier domain (Mialdun and Shevtsova, 2011b,a). The phase shift is then obtained by applying the inverse tangent function to the ratio of the imaginary part (Imₚ) and the real part (Reₚ), given by:

$$ \phi = \arctan \left( \frac{Im_p}{Re_p} \right) \qquad (1) $$

where φ is the phase shift. After application of the Fourier transform, the phase is initially obtained up to an arbitrary sign, which is inherent to the phase extraction procedure. We adjusted the sign by known distribution of refractive index.

At the beginning of the experiment, an interferogram of the initial state is recorded after the first step of cell filling. This interferogram
<figure><img src="image_3.png" /><figcaption>Fig. 2. Schematic illustration of the post-processing steps of an image in an experiment involving the superposition of 3.026 mol/L NaCl and 1.7 mol/L Na₂SO₄. a) Recorded interferograms (full field of view). b) Wrapped phase maps. c) Wrapped phase shift map, see Eq. (2). d) Refractive index map, see Eq. (4). Color bar depicts refractive index variations in the system.</figcaption></figure>

serves as the reference. Both the reference and subsequent interferograms are processed using the Fourier-transform method to obtain the their phase shift $\phi_r$. The phase shift difference $\Delta\phi$ is then computed as:

$$ \Delta\phi = \phi - \phi_r \qquad (2) $$

This subtraction also removes optical distortions inherent in the physical setup. The wrapped phase map $\Delta\phi$ is then processed to eliminate the 2$\pi$ phase ambiguity (branch cut discontinuity) inherent in the arctangent function, yielding the unwrapped phase difference:

$$ \Delta\varphi = \Delta\phi \pm 2\pi k \qquad (3) $$

where $k$ is an integer that accounts for the phase shift unwrapping ambiguity that results from the phase periodicity in the interferograms and $\Delta\varphi$ is the unwrapped optical phase shift difference. The phase unwrapping was performed using a two-dimensional phase-unwrapping algorithm implemented in Python, based on the algorithm developed by Herráez et al. (2002). The unwrapped phase shift difference is linked to the refractive index variation field $\Delta n$ in the diffusion cell by:

$$ \Delta n = \frac{\lambda}{2\pi l} \Delta\varphi \qquad (4) $$

where $\lambda$ is the laser wavelength, and $l$ is the optical path length of the cell. See Fig. 2 to view all the post-processing steps. The topology of the refractive index maps, as in panels (d), reveals the spatio-temporal evolution of the double diffusion instability in the system.

In the experiments, the interferometric signal is linked to variations in the refractive index of the solution. For incompressible liquids, refractive index variations can, in principle, originate from changes in pressure, temperature, and solute concentration. The experiments are conducted at ambient pressure, the room temperature is well controlled, and the two salts are not reactive. Thermal effects related to the mixing process are also assumed to be negligible. We therefore attribute the refractive index variations entirely to the evolving concentrations of NaCl and Na₂SO₄. The variation in refractive index is therefore related to variations in the concentration of the two salts by the relationship:

$$ \Delta n = \frac{\partial n}{\partial c_1} \Delta c_1 + \frac{\partial n}{\partial c_2} \Delta c_2 \qquad (5) $$

In practice, it is challenging to separate the contributions of both solutes in a purely optical measurement. A single-wavelength interferometer

provides access only to the combined refractive index field, and cannot unambiguously reconstruct the individual concentration fields in a ternary mixture. Even two-wavelength optical approaches suffer from severe ill-conditioning of the contrast-factor matrix, which makes determining individual components extremely inaccurate (Shevtsova et al., 2011). For this reason, we treat the refractive index as an integral quantity that reflects the combined effect of both concentrations, weighted by their contrast factors. Because of this experimental limitation, we cannot present our experimental results on instability in terms of the individual Rayleigh numbers of the solutes or their ratio, as is sometimes done in other studies (Kumar et al., 2018; Huppert and Moore, 1976; Huppert and Turner, 1981).

From a physical point of view, variations in refractive index are closely linked to density variations in aqueous salt solutions. For liquid mixtures, this link is formalized through the Lorentz-Lorenz relation. Because, in our experiments, density variations arise primarily from solutal concentration gradients while temperature is kept constant, the refractive index serves as a robust proxy for density. Consequently, interferometry provides a direct means to visualize and track the spatio-temporal evolution of density-driven instabilities. Although the method does not yield absolute concentration fields, it enables us to identify the onset, symmetry, and dynamics of the double-diffusive convection patterns, offering essential insight into the underlying instability mechanisms.

The accuracy of our interferometer can be assessed from the noise level of the phase maps, which is approximately 0.1 rad (≈ 2$\pi$/60). This value is consistent with the typical precision of interferometric measurements. When converted into refractive-index accuracy under our experimental conditions, it corresponds to ≈ 1 · 10⁻⁵. Such small fluctuations are far below the range of refractive-index variations relevant for the present instability dynamics and therefore are not observable in the figures.

## 2.3. Determination of salt diffusion coefficients

Determining reliable multicomponent diffusion coefficients for the NaCl + Na₂SO₄ + H₂O system forms a major scientific component of this work. In aqueous electrolyte systems, charge-induced ion-ion interactions lead to the electrostatic coupling of ionic fluxes that results in a coordinated movement of charged species to maintain local charge bal-
<caption>**Table 1**</caption>

The initial concentration of samples c in mol/L, the concentration ratio R = C<sub>bot</sub>/C<sub>top</sub>, the density ratio R<sub>p</sub> = ρ<sub>bot</sub>/ρ<sub>top</sub>, the buoyancy ratio R<sub>buo</sub> = β<sub>bot</sub>c<sub>bot</sub>/β<sub>top</sub>c<sub>top</sub>, the diffusivity ratio of the initial samples δ = D<sub>bot</sub>/D<sub>top</sub>, the main D<sub>kk</sub>/10<sup>-9</sup>m<sup>2</sup>/s and cross-diffusion D<sub>kl</sub>/10<sup>-9</sup>m<sup>2</sup>/s coefficients of the ternary mixture NaCl/Na<sub>2</sub>SO<sub>4</sub>/H<sub>2</sub>O used in the simulations. Horizontal dashed lines separate groups with equal concentration ratios. The remaining physical properties of the salts can be found in Supplementary materials.

<table><thead><tr><th></th><th>c<sub>NaCl</sub></th><th>c<sub>Na<sub>2</sub>SO<sub>4</sub></sub></th><th>R</th><th>R<sub>p</sub></th><th>R<sub>buo</sub></th><th>δ</th><th>D<sub>11</sub></th><th>D<sub>12</sub></th><th>D<sub>21</sub></th><th>D<sub>22</sub></th></tr></thead><tbody><tr><td>#1</td><td>3.026</td><td>1.700</td><td>0.56</td><td>1.071</td><td>1.465</td><td>0.296</td><td>1.01</td><td>0.38</td><td>0.052</td><td>0.485</td></tr><tr><td>#2</td><td>2.637</td><td>1.484</td><td>0.56</td><td>1.063</td><td>1.509</td><td>0.371</td><td>1.08</td><td>0.35</td><td>0.041</td><td>0.52</td></tr><tr><td>#3</td><td>0.500</td><td>0.280</td><td>0.56</td><td>1.014</td><td>1.664</td><td>0.595</td><td>1.62</td><td>0.23</td><td>-0.067</td><td>0.84</td></tr><tr><td>#4</td><td>1.700</td><td>1.700</td><td>1.00</td><td>1.125</td><td>2.404</td><td>0.306</td><td>1.088</td><td>0.29</td><td>0.043</td><td>0.52</td></tr><tr><td>#5</td><td>0.500</td><td>0.500</td><td>1.00</td><td>1.040</td><td>2.841</td><td>0.542</td><td>1.56</td><td>0.22</td><td>-0.054</td><td>0.78</td></tr><tr><td>#6</td><td>0.280</td><td>0.280</td><td>1.00</td><td>1.022</td><td>2.920</td><td>0.591</td><td>1.60</td><td>0.16</td><td>-0.089</td><td>0.87</td></tr><tr><td>#7</td><td>4.130</td><td>0.700</td><td>5.90</td><td>1.067</td><td>20.186</td><td>2.139</td><td>1.38</td><td>0.35</td><td>0.07</td><td>0.63</td></tr><tr><td>#8</td><td>4.130</td><td>1.245</td><td>3.32</td><td>1.009</td><td>5.786</td><td>2.740</td><td>1.26</td><td>0.40</td><td>0.13</td><td>0.53</td></tr></tbody></table>

ance. The relevance of these Coulomb effects has been explored in details under diffusion-dominated and macroscopic advection-dominated conservative and reactive transport conditions. This has led to the development of multicomponent diffusion theories and the concept of multicomponent ionic dispersion (Appelo and Wersin, 2007; Muniruzzaman and Rolle, 2019). The starting point of all these theories is the Nernst-Planck equation for the flux of ionic species in solutions. In the absence of an external electric field, the null current condition applies (Tabrizinejadas et al., 2021). Due to the conservation of the neutrality of the electric charge, the species associate and the diffusive behavior resembles that of molecules. As a result, the experimental measurements of the diffusion coefficients of salts provide an effective or apparent value, equivalent to that of neutral molecules. It is therefore possible to draw a parallel between Fick's formalism of the salts-molecules flux with the ionic species expression flux, and to write the apparent diffusion coefficients as a function of the self-diffusion coefficients of species and the thermodynamic correction factors.

If we consider the two dissolved electrolytes as molecular species, using Fick's formalism the diffusive fluxes of the two compounds in the direction opposite to their concentration gradients can be written as follows:

$$ \vec{J}_1 = -D_{11} \vec{V} c_1 - D_{12} \vec{V} c_2 \quad (6) $$

$$ \vec{J}_2 = -D_{21} \vec{V} c_1 - D_{22} \vec{V} c_2 \quad (7) $$

The D<sub>ij</sub> represent what are called the main diffusion coefficients and the D<sub>ij</sub> the cross diffusion coefficients. In Supplementary material S1 we have detailed how ternaries diffusion coefficients of NaCl and Na<sub>2</sub>SO<sub>4</sub> have been obtained. The geochemical software PHREEQC was used to obtain the self-diffusion coefficients D<sub>aq,i</sub> of species i. The required ionic activity coefficients γ<sub>i</sub> were calculated using the Pitzer equations, associated with reliable sets of specific interactions parameters for the species. The thermodynamic correction factor terms evaluation need the derivative of the species activity coefficients with respect to the concentration. The procedure is also detailed in the Supplementary material S1. For validation, the procedure was applied the concentrations of the ternary mixture for which diffusion coefficient measurements have been made in the literature, demonstrating that the accuracy of the diffusion coefficient determination is within 10%.

The final values of the main and cross diffusion coefficients for the electrolytes and corresponding to the equilibrium concentrations of the eight experimental samples are summarized in Table 1. These values were used as fixed input parameters for all nonlinear numerical simulations presented in following sections.

### 3. Results

We have studied, both experimentally and numerically, the instability in a vertical two-layer system composed of aqueous NaCl-H<sub>2</sub>O and Na<sub>2</sub>SO<sub>4</sub>-H<sub>2</sub>O solutions. The initial compositions of the aqueous solutions, the important parameters and the diffusion coefficients obtained in the previous section for all considered cases are listed in Table 1. Here we present eight experiments, which are grouped as follows:

1. The molar concentration ratio R = C<sub>b</sub>/C<sub>t</sub> = 0.56 is maintained R = 0.56, while the concentration of each solution is varied.

2. Equal molar concentrations in both layers, R = 1, while the concentration of each solution is varied.

3. Salt positions inverted (NaCl-H<sub>2</sub>O at the bottom, Na<sub>2</sub>SO<sub>4</sub>-H<sub>2</sub>O at the top) with R = 3.3 and R = 5.9 to maintain gravitational stability.

In the first two groups, Na<sub>2</sub>SO<sub>4</sub>-H<sub>2</sub>O is on top and Na<sub>2</sub>SO<sub>4</sub>-H<sub>2</sub>O is on the bottom. In the third group, the location of the solutions is reversed: Na<sub>2</sub>SO<sub>4</sub>-H<sub>2</sub>O is on the top and NaCl-H<sub>2</sub>O is on bottom. The compositions are selected to keep the system gravitationally stable. Note, the concentration ratio R in the third group is much higher than in previous groups, i.e. 3.3 and 5.9.

Within each group, cases are ordered by decreasing density ratio R<sub>p</sub>. In contrast, the diffusivity ratio δ is increasing in this order. In Groups 1 and 2, R<sub>buo</sub> > 1 and δ < 1, hence a double-diffusive layer convection (DLC) instability is expected. In Group 3, δ > 1, leading to different buoyancy-driven dynamics.

#### 3.1. Experimental observations

##### 3.1.1. Asymmetrical and delayed onset of diffusive-layer-convection (DLC) instability

We begin by discussing the first group of experiments. In all three cases, the double-diffusive layer convection (DLC) instability was observed in a systems composed of NaCl-H<sub>2</sub>O and Na<sub>2</sub>SO<sub>4</sub>-H<sub>2</sub>O. The novelty of our findings lies in the experimental evidence of both a "delayed" and "asymmetrical" onset of instability.

For the first time, we observed an asymmetrical and delayed onset of double-diffusive convection (DLC) in a two-layer NaCl-Na<sub>2</sub>SO<sub>4</sub> system. This instability initiates in the upper layer and only later develops in the lower layer - a sequence that has not been reported in previous works. The delay in onset, which varies systematically with the diffusivity ratio δ, highlights the interplay between solute transport rates and interface sharpening. Interestingly, a possibility of delayed instability was anticipated theoretically by Trevelyan et al. (2011) in their linear stability analysis for δ < 1 < R<sub>buo</sub>. Our findings provide experimental confirmation of this prediction, while also revealing new asymmetries and plume dynamics that extend beyond the scope of prior analysis.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Experimentally observed DLC in Group 1, shown in terms of the refractive index (RI) measured by interferometry, which serves as a surrogate for the density field through the Lorentz-Lorenz relation. The first column shows the RI profile at the same times as the following snapshots (1, 2, 3 and 4, respectively). The subsequent snapshots illustrate the initial convection in the upper layer, followed by the development of bottom-driven convection. The non-linear color-maps were employed for enhancing the visibility of plumes in continuous 2-D refractive index distributions. To smooth out the local variations caused by instability, 1-D profiles (first panels) were obtained by averaging 2-D maps along horizontal direction.</figcaption></figure>

<caption>**Table 2**</caption>

Onset times of the DLC instability in the upper and lower layers for Groups 1 and 2, obtained from interferometric measurements. Times are given in seconds, with an uncertainty of ±10 s.

<table><thead><tr><th rowspan="2"></th><th colspan="3">Group 1</th><th colspan="3">Group 2</th></tr><tr><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th></tr></thead><tbody><tr><th>Top</th><td>100</td><td>0</td><td>250</td><td>160</td><td>730</td><td>∞</td></tr><tr><th>Bottom</th><td>170</td><td>350</td><td>400</td><td>700</td><td>∞</td><td>∞</td></tr></tbody></table>

Quantitative data on the onset time is given in Table 2. The onset time is defined as the first instant at which a persistent deformation of the initially flat diffusion interface is observed, marking the transition from pure diffusion to sustained convective motion. This was determined by frame-by-frame analysis of the video recordings using the same criterion for both layers. The uncertainty of ±10 s reflects the temporal resolution of the recordings and the gradual nature of the transition.

The term “asymmetrical onset” refers to convection first appearing in the upper layer while the bottom layer remains quiescent. This behavior is due to the faster diffusion of NaCl ($D_{11}$, upper layer) compared to $Na_2SO_4$ ($D_{22}$, lower layer). Consequently, NaCl diffuses downward quickly, creating a localized density inversion right below the interface, whereas $Na_2SO_4$ diffuses upward much more slowly. The development of convection over time for all six cases in Groups 1 and 2 is shown in Fig. 3, expressed in terms of the refractive index.

As NaCl diffuses out of the region just above the interface (into the bottom layer), it leaves behind a region with reduced NaCl concentra-

tion. Since $Na_2SO_4$ does not diffuse upward efficiently, this zone does not get compensated by solute from below. This creates a locally low density region just above the interface which is shown by blue color. This blue zone is not yet convective, it's diffusion-driven depletion. The density inversion between this depleted region (light) and the fluid above it (denser NaCl) causes DLC.

By contrast, the slow diffusion of $Na_2SO_4$ into the NaCl layer produces a more gradual modification of the density profile in the lower part of the diffusion zone. In order to become unstable, the sublayer below must reach a certain critical thickness and density contrast threshold. The local density gradient at the bottom only becomes unstable when enough NaCl has diffused downward. Later on, the flow in the bottom part becomes dominant and overtakes the initial top-driven DLC. The refractive index profile at $t=1000$ s, shown by curve 4 in the first panel of Fig. 3 (indirectly indicating density) highlights the strong convection in the bottom layer at later times.

Taking the time delay of the instability into account, we can estimate the distance at which the depleted layer forms. Given that $\tau \sim l^2/D_{eff}$, the distance is $l \sim \sqrt{\tau D_{eff}}$. The diffusion coefficient of NaCl in water is approximately $1.5 \times 10^{-9} \text{ m}^2/\text{s}$ at 25 °C. Therefore, in the Group 1, for a delay of 100 s, the depletion zone will be located at $l \sim 0.5 \text{ mm}$ above interface and for a delay of 150 s, $l \sim 1 - 1.5 \text{ mm}$ below. Actually, the experimental observations in Fig. 3 confirm these estimations.

A further insight comes from the diffusivity ratio $\delta = D_{bot}/D_{top}$, which determines the extent to which the gradient is smoothed by differential diffusion. When $D_{bot} \ll D_{top}$, the top-layer solute diffuses rapidly downward, effectively smearing out the compositional interface before any instability can grow. This smearing delays the onset of convection.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Experimentally observed DLC in Group 2 (cases 4, 5 and 6), shown in terms of the phase difference from interferometric measurements. Cases #4 and # 5 illustrated delayed asymmetrical instability, whereas case # 6 (bottom row) did not reveal any instability. Wrapped phase visualization enhances sensitivity to small refractive index variations, revealing subtle convective patterns in Case #4. Remind that the field of view is 11 mm × 11 mm.</figcaption></figure>

By the time plumes develop, the concentration gradient has become more diffuse. This results in fewer, thicker plumes with a reduced vertical extension. As Table 1 shows, δ increases from case #1 (0.296) to case #3 (0.595), while R<sub>buo</sub> changes only modestly (1.465 → 1.664). Fig. 3 demonstrates that case #2, with an intermediate δ and R<sub>buo</sub>, exhibits the most vigorous convection: plumes appear almost immediately, grow tall, merge, and then coarsen into a few wide cells by 1000 s. Curve 3 (green, t=355s) in the first panel of Fig. 3 underlines the significant concentration gradient just above the initial contact line. In contrast, case #3 (highest δ) shows the most delayed onset (first fingers at ≈ 250 s in Table 2) and only a few thick, shallow plumes, consistent with the over-smearing effect.

Similarly, case #1 (lowest δ) shows a weak and delayed convection, attributed to a slow gradient build-up despite its large buoyancy forcing. Together, Table 2 and Fig. 3 confirm that an intermediate diffusivity ratio, as in case #2, offers the optimal balance between gradient sharpness and buoyant forcing, maximizing DLC activity.

In the second group of experiments, the buoyancy ratio $R_{\text{buo}}$ is approximately twice as large as in the first group, ranging from 2.40 to 2.84. At first glance, this stronger buoyancy would be expected to promote earlier and more vigorous instability. However, the experiments reveal a different behavior: DLC occurs only in Case #4, and with a significant delay in case #5 (~160 s at the top and ~730 s at the bottom), while no convection is observed in Case #6. The evolution of diffusive patterns in Group 2 is shown in Fig. 4 in terms of the wrapped phase for better visibility of the subtle patterns, particularly in case #5.

### 3.1.2. Double-diffusive (DD) instability

In the third group of experiments, the positions of the salts are reversed compared to previous configurations: NaCl-H₂O is placed at the bottom and Na₂SO₄-H₂O at the top. In this configuration, the lower component, NaCl, diffuses faster than the upper one, leading to the development of a double-diffusive (DD) instability, typically manifesting in the form of salt fingers. Both systems are gravitationally stable though the density ratios are small, $R_ρ = 1.067$ and $R_ρ = 1.009$, respectively. Despite the similar arrangement, Cases #7 and #8 show markedly distinct behaviors.

To present the results of case #7, we use in Fig. 5 a wrapped phase map instead of a refractive index map. If refractive index maps were used instead, the overall structure would be less distinguishable, although minor undulations would still be visible near the bottom layer. The wrapped phase maps reveal a long diffusion-controlled regime, during which minor undulations gradually develop at the leading edge.

These undulations appear as a fine waviness that becomes more pronounced with time, and may give the visual impression of a staircase-like structure in the phase images. However, this apparent layering is not reflected in the actual refractive index profiles, which remain smooth over the entire domain as shown in the first panel. The visual effect arises from phase wrapping of the optical signal, not from the formation of well-defined convective layers.

Unlike case #7, the results of case #8 are presented in the form of a refractive index map to make the details easier to see. Black/white are used here for to avoid conflicting views between cases, although a color map was used for the refractive index in Fig. 3. Case #8 exhibits an immediate onset of fingering upon contact between the layers. These fingers are unusually fine with minimal horizontal extent and low amplitude. This makes them appear more like a faint, vertically striated texture than distinct convective fingers. The growth rate is significantly slower than those seen in DLC scenarios. This weak behavior can be attributed to the extremely small density ratio ($R_ρ = 1.009$), which places the system very close to neutral stability. Additionally, the large diffusivity ratio between the two mixtures ($δ = D_{\text{NaCl}}/D_{\text{Na}_2\text{SO}_4} = 2.74$) implies a strong asymmetry in mass transport: the lower, faster-diffusing component promotes interface destabilization, while the upper component lags behind in response. This imbalance may contribute to the suppression of coherent convective structures. At later times (e.g., t=950 s), the structure becomes increasingly irregular, lacking the periodic layering observed in case #7. These fingers are likely salt-rich plumes of NaCl penetrating upward, destabilizing the fluid above, and being modulated by slow diffusion of Na₂SO₄ downward. The dynamics appear to enter a quasi-chaotic regime, characterized by unsteady, spatially incoherent convective motions, where interactions between fingers lead to lateral merging, detachment, and chaotic evolution. This is a hallmark of nonlinear saturation of diffusion layer convection in systems close to gravitational neutrality.

Somewhat similar patterns have been observed in numerical simulations (Radko, 2005) and experiments (Huppert and Turner, 1981; Turner and Stommel, 1964), particularly in systems with small density ratios and strong asymmetry in solute diffusivities. These studies show that when stabilizing and destabilizing effects are nearly balanced, the system can evolve into a quasi-chaotic state characterized by intermittent layering, coalescing fingers, and irregular flow, as observed in our refractive index maps for case #8.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Evolution of convection in the third group of experiments shown using wrapped phase images. Panels in upper row correspond to Case #7, showing the development of salt fingers and the eventual formation of a staircase pattern after t = 600 s. Wrapped phase visualization fits better to case # 7. Panels in the lower row show case # 8, using refractive index visualization to enhance sensitivity and reveal subtle convective patterns. Fine-scale vertical striations appear immediately upon contact, with a staircase-like structure emerging after t = 950 s.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Results of the nonlinear simulations for group #1 showing the evolution of the NaCl inside the full cell (17 mm × 58 mm) correspond to two time points: t=175 s, during active plume growth, and t=350 s, when the leading plume approaches the upper boundary. The comparison highlights differences in the onset timing and plume morphology across Cases 1-3. The diffusion coefficients are taken from Table 1.</figcaption></figure>

## 3.2. Nonlinear simulations

In previous sections, we explained how diffusion coefficients are determined using a combination of thermodynamic estimates and numerical calculations. Using a Mach-Zehnder interferometer, we conducted experiments and observed the onset and evolution of convective regimes. In this section, we present nonlinear simulations that use these diffusion coefficients as input to reproduce the experimental observations. These simulations are not merely illustrative; they act as a stringent test of our overall approach. By accurately reproducing the key features seen in the experiments, it will be confirmed that the transport properties were correctly identified.

At the macroscopic scale, the spatio-temporal dynamics of the diffusive system is described by the Navier-Stokes and mass transfer equations in the Boussinesq approximation. The mass transfer equations, written for each solution, include diffusion and cross-diffusion coefficients. The Boussinesq approximation assumes that the density fluctuations that occur with the onset of fluid motion are caused solely by concentration effects. All variations in fluid properties other than density are neglected.

$$ \mathrm{div} \, \vec{V} = 0, \qquad (8) $$

$$ \frac{\partial \vec{V}}{\partial t} + \vec{V} \cdot \nabla \vec{V} = - \frac{1}{\rho_r} \nabla p + \nu \Delta \vec{V} + \vec{g} \rho^* \qquad (9) $$

$$ \frac{\partial w_1}{\partial t} + \vec{V} \cdot \nabla w_1 = D_{11} \Delta w_1 + D_{12} \Delta w_2 \qquad (10) $$
<figure><img src="image_8.png" /><figcaption>Fig. 7. (The results are only displayed in the experimental field of view (11 mm × 11 mm).) Results of the nonlinear simulations for Group 2 showing the evolution of the NaCl corresponding to two time points: at t = 350 s and at t = 700 s.</figcaption></figure>

$$ \frac{\partial w_2}{\partial t} + \vec{V} \cdot \nabla w_2 = D_{21} \Delta w_1 + D_{22} \Delta w_2 \quad (11) $$

Here $w_i$ are the mass fractions, and $\vec{V}$ is the velocity; $p$ is the difference between total and hydrodynamic pressure; $v$ is the kinematic viscosity and $\rho^* = \rho - \rho_r$.

$$ \rho = \rho_r [1 + \beta_{w1}(w_1 - w_{1r}) + \beta_{w2}(w_2 - w_{2r})], \quad (12) $$

$$ w_{ir} = (w_{i0}^t + w_{i0}^b)/2 \quad \beta_{c_i} = \frac{1}{\rho} \frac{\partial \rho}{\partial w_i} \Big|_{c_{i0}} $$

where $\rho_r = (\rho_0^t + \rho_0^b)/2$ are the reference values. Information on the physical properties of the mixtures, as well as the correspondence between molar concentration and mass fraction units, are provided in the supplementary materials.

Hereafter, the notations $D_{ij}$ and $D_{il}$ will be used for the full set and the main terms of the diffusion matrix, respectively. The values of the diffusion coefficients are taken from Table 1.

Heat and mass transfer equations are typically formulated in the mass reference frame, and we follow this convention here, although the diffusion coefficients used were obtained in the solvent reference frame. Previous studies (Kozlova et al., 2019) have shown that for mixtures behaving nearly ideally, the choice of reference frame has little impact on the values of $D_{ij}$. Given that the diffusion coefficients were determined with an uncertainty of about 10%, this margin covers the expected differences between the mass and solvent frames of reference.

Eqs. (8)–(12) were solved within the OpenFOAM framework, enabling detailed numerical simulations of instability mechanisms in multicomponent systems. The numerical code has been validated in earlier studies (Šeta et al., 2019; Šeta, 2020).

Below we present evolution of convection in terms of the mass fraction field $w_1$ (NaCl). The evolution of the full set of variables can be seen in Supplementary Movie 1, where we selected Case 2 as the representative example. Fig. 6 shows the evolution of $c_1$ (NaCl) in Group 1 at two selected time points: $t = 175$ s, during the active growth phase of the instability, and $t = 350$ s, when the most advanced plume reaches the upper boundary.

The comparison between the cases reveals that the instability develops more rapidly in Case 2 than in Case 1, even though their density ratios $R_\rho$ are nearly identical. This trend aligns well with the experimental observations. A more detailed agreement is observed in the morphology of the plumes: Case 2 displays the thinnest plumes, followed by moderately thicker plumes in Case 1, and significantly broader structures in Case 3. This variation correlates with the diffusivity ratio $\delta$; as $\delta$ increases, the plumes become wider.

Fig. 7 presents the simulated evolution of the NaCl mass fraction field $w_1$ for Group 2 at two selected time points: $t = 350$ s, when the first convective plumes begin to form, and $t = 700$ s, when the leading plume reaches the upper edge of the experimental field of view.

In this group, the initial concentrations at the top and bottom boundaries are equal, resulting in a markedly weaker and more delayed development of instability compared to Group 1. Notably, in Case 5, the onset of convection is significantly delayed, while in Case 6, no instability

develops within the observed time window. Consequently, the presentation of the numerical results was restricted to match the experimental field of view, focusing on direct comparability.

The first round of simulations used diffusion coefficients from Table 1. These initial results produced mixed agreement with the experiments. Case 5 aligned well with the experimental observations, showing a delayed but eventual onset of convection. However, Case 4 remained entirely stable in the simulation, despite the presence of convection in the experiment. Conversely, Case 6 developed instability in the simulation, while it remained stable in the experimental setup.

To reconcile these discrepancies in Cases 4 and 6, we slightly adjusted the diffusion coefficients within the bounds of the estimated 10% measurement uncertainty. In Case 4, increasing $D_{11}$ as (1.088 → 1.15) × 10⁻⁹m²/s and decreasing $D_{22}$ as (0.52 → 0.47) × 10⁻⁹m²/s shifted the simulated system from a stable to an unstable regime, matching experimental behavior. In Case 6, adjusting $D_{11}$ as (1.6 → 1.4) × 10⁻⁹m²/s and $D_{22}$ as (0.87 → 1.0) × 10⁻⁹m²/s had the opposite effect-stabilizing the system in line with the experimental result. In the geometrical representation, these modifications can be seen later in Fig. 9.

The corrections, though minor, are fully justified given the inherent uncertainty in diffusion coefficient estimation. The improved simulations accurately capture both the onset of the instability and the structural characteristics of the convective patterns. This underscores the sensitivity of the system to transport parameters and reinforces the validity of the experimental and numerical framework used in this study.

Interestingly, Bratsun et al. (2022) reported noticeably different stability behaviors in two aqueous systems when the concentration ratio was kept constant but the individual concentrations were varied. In the NaNO₃-H₂SO₄ system, the solution remained stable at low concentrations but became unstable upon a proportional increase in concentration. Conversely, for LiCl and NaNO₃ solutions, higher concentrations led to enhanced stability. In contrast, our system exhibits a non-monotonic relationship between concentration and stability, as demonstrated in the preceding analysis. This suggests that the dependence of convective behavior on composition is more complex and system-specific.

The results of the non-linear simulations of Group 3 are fully consistent with the experimental observations, as can be seen in Fig. 8. In case #7, the system undergoes diffusion over a long period of time. The first clearly visible fingers appear at $t \sim 500$ s (not shown) and evolve slowly. For this reason, the second snapshot in Fig. 8 presents the pattern at a later time, $t = 1000$ s. Over the next 1000s, the fine, hair-like structures made little progress in their evolution.

In case 8, the simulated fingers resemble narrow, wavering filaments, consistent with experimental observations. This creates the impression of a slow, vertical shimmering rather than discrete convection patterns. The development of this pattern is notably slow; even at 1750 seconds, it remains confined to less than half the vertical extent of the cell. Another notable feature is the continuous generation of fine fingers at the interface, which gradually mix and merge some distance above it, forming a nearly uniform concentration field. This is evident in the snapshot
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Results of the nonlinear simulations for Group 3, showing the evolution of NaCl within the entire experimental cell (17 mm × 58 mm) at representative time points for two cases: case 7 and case 8.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Stability maps for the first six cases. The upper row shows maps for Group 1 with the concentration ratio $R=0.56$ while the lower one for Group 2 with $R=1$. The crosses on the maps indicate the initial state of the system. For both groups the density and buoyancy ratio is decreased from the left to right. In the space below the red diagonal line, the system is prone to Rayleigh-Taylor instability; the blue and yellow regions correspond to DLC and DD, respectively; the rest of the parameter space corresponds to the stable state. The dashed lines in cases #4 and #6 show how diffusion coefficients were adjusted in non-linear simulations. Each panel corresponds to a specific set of diffusion coefficients from Table 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

at $t=1750$ s, where the base appears distinctly darker. This behavior is consistent with the expected dynamics of double-diffusive convection in systems close to neutral stability, where the stabilizing density gradient is too weak ($R_p = 1.009$) to suppress differential diffusion effects (Radko, 2005).

These results indicate that the simulations successfully reproduce not only the onset timing of the instability but also the structural differences in convective patterns. This confirms the reliability of the estimated

transport coefficients and the effectiveness of the numerical model in capturing the key physical mechanisms observed experimentally.

### 3.3. Additional criterion of instability in ternary systems

To gain further insight into the distinct behavior exhibited by the different groups, we now turn to the analytical solution of the 1D diffusion problem in a ternary system with constant diffusion coefficients. The diffusion equations in three-components system are written
as

$$ \frac{\partial w_1}{\partial t} = D_{11}\Delta w_1 + D_{12}\Delta w_2, \qquad (13) $$

$$ \frac{\partial w_2}{\partial t} = D_{21}\Delta w_1 + D_{22}\Delta w_2. \qquad (14) $$

where the mass fractions $w_j$ are functions of position $z$ and time $t$, and $D_{ij}$ are the diffusion coefficients. At the initial time of contact ($t=0$), the two solutions are separated by a sharp interface at $z=0$, and the density profile follows a step function. The initial composition differences between the layers are denoted as $(\Delta w_{01}, \Delta w_{02})$.

For a purely diffusive process with constant diffusion coefficients, Eqs. (13) and (14) can be solved analytically (Šeta et al., 2019; Šeta, 2020). These solutions delineate stability regions corresponding to different system behaviors depending on the initial concentration differences: unstable, stable, double-diffusive instability (DD) in the form of fingers, and diffusion-layer convection (DLC). This analytical framework fully characterizes the system at its initial stage and serves as a predictor of its subsequent evolution.

The solutions can be represented in terms of a stability map plotted in $(\Delta w_{01}, \Delta w_{02})$ space, which is equivalent to $(\Delta \rho_{01}, \Delta \rho_{02})$ space. In this map, Rayleigh-Taylor-type instability occurs when $R_{buo} < 1$, i.e., when $\rho_{top} > \rho_{bottom}$, and this region lies below the red diagonal line. The ratio $\beta_{c1}\Delta c_1/\beta_{c2}\Delta c_2$ is mathematically equal to the ratio of the solutal Rayleigh numbers $Ra_{c1}/Ra_{c2}$, and the regimes could also be expressed directly in Rayleigh-number space. We present the results in terms of this ratio because it is the conventional parameter used in the double-diffusive literature (Bratsun et al., 2022; Mcdougall, 1983), which ensures consistency with prior work and facilitates comparison with existing studies.

The stability maps for Groups 1 and 2, corresponding to the first six experimental cases, are shown in Fig. 9. In these maps, the blue region corresponds to DLC instability and yellow to DD instability. The extent of the DLC region (shaded blue) increases as the diffusivity ratio $\delta$ rises. Based on the positions of the crosses on the maps, diffusion-layer convection (DLC) is predicted for all cases from #1 to #6. Our experiments reveal that when the buoyancy ratio $R_{buo}$ varies only slightly between cases within the same group, the strength of the instability is influenced not only by the diffusivity ratio $\delta$, but also by the initial position of the system on the stability map. When the system is close to the boundaries of the DLC region, it becomes more sensitive to perturbation growth due to reduced stability margins.

This trend is clearly seen within Group 1 in Fig. 9. For example, instability evolves more quickly in Case #2 than in case #1, even though their density and buoyancy ratios are nearly identical. This suggests that the diffusivity ratio $\delta$ accelerates instability growth. However, the trend is not purely monotonic: case #3, which has a $\delta$ twice as large as case #1, does not exhibit faster evolution compared to case #2.

This apparent discrepancy can be explained by the initial positioning on the stability map. In case #2, the cross is located deeper within the DLC region, where perturbations grow faster. By contrast, cases #1 and #3 lie closer to the DLC region boundary, where the system remains marginally stable for a longer time before instability sets in. Thus, for a given $R_{buo}$, the development of instability is governed by both the diffusivity ratio, $\delta$, and the initial position relative to the stability boundaries.

Observations from Group 2 further support this interpretation. Although the buoyancy ratio $R_{buo}$ in Group 2 is approximately twice that of Group 1, and the corresponding diffusivity ratios are comparable or even higher, the onset of DLC is noticeably delayed (see Table 2). This suggests that the instability threshold is not fully governed by local buoyancy forcing and diffusivity contrast alone. Rather, it seems that the critical factor is the system's initial position within the stability map. This determines how close it is to the boundaries of instability and thus how long it takes to trigger convective motion.

In case #4, the cross is located directly on the DLC boundary, bordering a stable region. Instability does develop but only after a moderate

delay. In case #6, the cross lies near the corner of the map, close to the boundary between stable, DLC, and DD regions. Here, instability is absent within the observation window. Case #5, with a slightly more favorable initial position, develops instability but still with a substantial delay compared to Group 1.

Fig. 10 shows that all cases in Group 3 fall within the yellow region on the stability diagram, indicating susceptibility to double-diffusive fingering instability. In case #7, the initial conditions lie deep within the fingering regime. However, the experimental results show a delayed onset of fingering, occurring only at the leading edge of the concentration front. This delayed behavior is attributable to the small density ratio and a high buoyancy ratio (approximately 20 times larger at the bottom), which act as stabilizing factors suppressing early fingering activity.

In case #8, the initial condition is located close to the neutral stability curve. Systems near this threshold are only marginally unstable. Hence, the observed fingering behavior, characterized by low velocity, fine spatial scales, and minimal contrast, is expected and consistent with weak solutal convection close to marginal stability.

Our results show that the initial compositions play a key role in determining the position of the system on the stability map (i.e., cross). This position ultimately determines the onset and nature of instability, highlighting the importance of composition tuning in predicting and controlling double-diffusive convection.

## 4. Conclusions

We investigated double-diffusive (DD) and diffusion-layer convection (DLC) in two characteristic salts of saline aquifers, i.e., superimposed aqueous solutions of NaCl and Na₂SO₄. This study combined the determination of diffusion coefficients in ternary mixtures, experimental observations using a Mach-Zehnder interferometer, analytical interpretation, and nonlinear numerical simulations.

To determine the diffusion coefficients of the ternary mixture, we treated the two dissolved electrolytes as molecular species within the framework of Fick's formalism. The geochemical software PHREEQC, which solves the full Nernst-Planck equation numerically, was used to obtain the self-diffusion coefficients of all ionic species. The required ionic activity coefficients for computing the main and cross-diffusion coefficients, were calculated using the Pitzer equations, associated with reliable sets of specific interactions parameters for the species. For validation, the procedure was applied to a well-characterized reference system, demonstrating that the accuracy of the diffusion coefficient determination is within 10%.

Optical digital interferometry based on a Mach-Zehnder interferometer was employed to track minute concentration changes in a 1 mm thick diffusion cell. The recorded fringe patterns were processed to reconstruct the spatial distribution of the optical phase. The resulting data were interpreted through either wrapped phase maps or refractive index fields. Non-linear color maps were applied to the reconstructed 2-D refractive index distributions to enhance the visibility of plumes, while additional post-processing techniques improved the accuracy of the concentration profiles.

The diffusion coefficients obtained from the ternary analysis were used as input for nonlinear simulations, which successfully reproduced the observed experimental dynamics. This agreement validates both the experimental methodology and the theoretical framework, providing a coherent and self-consistent picture of the underlying physical mechanisms.

To ensure a consistent scientific outcome in a system governed by many interdependent parameters, we categorized the scenarios into three groups, each containing two or three cases. Within each group, the concentration ratio between NaCl and Na₂SO₄ was held constant, while the absolute concentrations were varied.

In the first two groups, where the faster-diffusing NaCl was placed in the upper layer, we observed diffusion-layer convection (DLC) when instability occurred. A key novelty of our findings is the experimental
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Stability maps for the third group. In the space below the red diagonal line, the system is prone to Rayleigh-Taylor instability; the blue and yellow regions correspond to DLC and DD, respectively; the rest of the parameter space corresponds to the stable state. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

demonstration of both a delayed and asymmetric onset of convective instability.

In the third group, the positions of the two salts were reversed, resulting in double-diffusive instability manifested as vertically aligned fingers. Their growth was significantly slower than in the DLC scenarios, likely due to a very small density ratio. Both the experiments and simulations revealed an unusual finger morphology: the structures are exceptionally fine and resemble a faint, vertically striped texture rather than distinct convective fingers. Notably, they do not disrupt the interface but are instead continuously generated there, gradually mixing and merging some distance above and below the initial boundary to produce a nearly uniform concentration field.

The studied cases spanned a wide range of density and diffusivity ratios, yet classical stability criteria did not fully capture the observed behavior. We found that an additional parameter - the system's initial position on the stability map - must be taken into account. This position is determined by the full diffusion matrix, underscoring the importance of multicomponent diffusion effects in accurately predicting the onset of instability in ternary systems.

Our results provide new insight into solute transport in saline aquifers, highlighting the need for refined stability analysis that accounts for cross-diffusion and thermodynamic interactions. Further work could focus on improving the accuracy of the calculated diffusion coefficients and their dependence on composition (Šeta et al., 2025), so that the simulations can better guide the experiments.

### Supporting information

The important part of the work, determination of the ternary diffusion coefficients of the mixture NaCl + Na₂SO₄ + H₂O is detailed in the first section of the supplementary materials. The second section presents correspondence between molar concentration and mass fraction units, information on the physical properties of the mixtures, as well as binary diffusion coefficients from various literature sources, are provided in the supplementary materials. The time evolution of NaCl concentration during the mixing of the solutions for the case 2, as a representative example, can be seen in Supplementary Video 1.

### CRediT authorship contribution statement

**Christian O. Oko:** Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis; **Berin Šeta:** Writing – review & editing, Visualization, Validation, Software, Investigation, Formal analysis, Conceptualization; **Aliaksandr Mialdun:** Writing – review & editing, Visualization, Investigation, Formal analysis, Data curation; **Mounir M. Bou-Ali:** Writing – review & editing, Investigation, Funding acquisition, Conceptualization; **Valentina Shevtsova:** Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Conceptualization; **Arnault Lassin:** Writing – review & editing,

Software, Methodology, Investigation, Formal analysis; **Henri Bataller:** Writing – original draft, Supervision, Software, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgement

The work is supported by Elkartek program (KK-2023/00041- MMASINT) and Research Group Program (IT1505-22) of the Basque Government, and Grant PID2021-124232OB-I00 (Treated, AEI/FEDER, UE), funded by MCIN/AEI/10.13039/501100011033 of the Spanish Government and, as appropriate, by "ERDF A way of making Europe", by the "European Union". C.O.O and H.B. would like to acknowledge the financial support from the French Centre National d'Etudes Spatiales (CNES), the E2S-UPPA Hub Newpores project, the Industrial Chair CO2ES supported by the Investissements d'Avenir French program managed by ANR (No. ANR16IDEX0002) and the Nigerian Petroleum Technology Development Fund (PTDF).

### Supplementary material

Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ces.2025.123224

### References

* Albright, J.G., Gillespie, S.M., Rard, J.A., Miller, D.G., 1998. Ternary solution mutual diffusion coefficients and densities of aqueous mixtures of nacl and na2SO4 at 298.15 k for six different solute fractions at a total molarity of 1.000 mol·dm⁻3. J. Chem. Eng. Data 43, 668–675. https://doi.org/10.1021/je9800196

* Annunziata, O., Rard, J.A., Albright, J.G., Paduano, L., Miller, D.G., 2000. Mutual diffusion coefficients and densities at 298.15 k of aqueous mixtures of nacl and na2SO4 for six different solute fractions at a total molarity of 1.500 mol·dm⁻3 and of aqueous na2SO4. J. Chem. Eng. Data 45, 936–945. https://doi.org/10.1021/je000134s

* Appelo, C.A.J., Wersin, P., 2007. Multicomponent diffusion modeling in clay systems with application to the diffusion of tritium, iodide, and sodium in opalinus clay. Environ. Sci. Technol. 41, 5002–5007. https://doi.org/10.1021/es0629256

* Applin, K.R., Lasaga, A.C., 1984. The determination of SO42-, naSO4-, and mgSO40 tracer diffusion coefficients and their application to diagenetic flux calculations. Geochim. Cosmochim. Acta 48, 2151–2162. https://doi.org/10.1016/0016-7037(84)90395-8

* Archer, D.G., 1992. Thermodynamic properties of the nacl + h2o system. II. thermodynamic properties of nacl(aq), nacl·2h2(cr), and phase equilibria. J. Phys. Chem. Ref. Data. 21, 793–829. https://doi.org/10.1063/1.555915

* Bratsun, D.A., Oschepkov, V.O., Mosheva, E.A., Siraev, R.R., 2022. The effect of concentration-dependent diffusion on double-diffusive instability. Phys. Fluids 34 (3), 034112. https://doi.org/10.1063/5.0079850

* Fu, J., Paduano, L., Rard, J.A., Albright, J.G., Miller, D.G., 2001. Mutual diffusion coefficients and densities at 298.15 k of aqueous mixtures of nacl and na2SO4 at high
concentrations with nacl solute fractions of 0.9000. J. Chem. Eng. Data 46, 601–608.
https://doi.org/https://doi.org/10.1021/je000357i

Fu, J., Rard, J.A., Paduano, L., Albright, J.G., Miller, D.G., 2002. Mutual diffusion coefficients and densities at 298.15 k of aqueous mixtures of nacl and na2SO4 with nacl solute fractions of 0.9500, trace diffusion coefficients of SO42- in nacl(aq), and trace refractive index increments and partial molar volumes of na2SO4 and nacl. J. Chem. Eng. Data 47, 496–512. https://doi.org/https://doi.org/10.1021/je010304g

Hall, J.R., Wishaw, B.F., Stokes, R.H., 1953. The diffusion coefficients of calcium chloride and ammonium chloride in concentrated aqueous solutions at 25°C. J. Am. Chem. Soc. 75, 1556–1560. https://doi.org/10.1021/ja01103a011

Herráez, M.A., Burton, D.R., Lalor, M.J., Gdeisat, M.A., 2002. Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path. Appl. Opt. 41 (35), 7437–7444. https://doi.org/10.1364/AO.41.007437

Huppert, H.E., Moore, D.R., 1976. Nonlinear double-diffusive convection. J. Fluid Mech. 78 (4), 821–854. https://doi.org/10.1017/S0022112076002759

Huppert, H.E., Turner, J.S., 1981. Double-diffusive convection. J. Fluid Mech. 106, 299–329.

Kozlova, S., Mialdun, A., Ryzhkov, I., Janzen, T., Vrabec, J., Shevtsova, V., 2019. Do ternary liquid mixtures exhibit negative main fick diffusion coefficients? Phys. Chem. Chem. Phys. 21, 2140–2152. https://doi.org/10.1039/C8CP06795C

Kumar, V., Srivastava, A., Karagadde, S., 2018. Compositional dependency of double-diffusive layers during binary alloy solidification: full-field measurements and quantification. Phys. Fluids 30, 113603. https://doi.org/10.1063/1.5049135

Luo, A., Li, Y., Chen, X., Zhu, Z., Peng, Y., 2022. Review of CO2 sequestration mechanism in saline aquifers. Natur. Gas Indus. B 9 (4), 383–393. https://doi.org/https://doi.org/10.1016/j.ngib.2022.07.002

Mcdougall, T.J., 1983. Double-diffusive convection caused by coupled molecular diffusion. J. Fluid Mech. 126, 379–397. https://doi.org/10.1017/S0022112083000221

Mialdun, A., Shevtsova, V., 2011a. Digital interferometry as a powerful tool to study the thermodiffusion effect. CR Mécanique 339 (5), 362–368.

Mialdun, A., Shevtsova, V., 2011b. Measurement of the soret and diffusion coefficients for benchmark binary mixtures by means of digital interferometry. J. Chem. Phys. 134 (4), 044524. https://doi.org/10.1063/1.3546036

Muniruzzaman, M., Rolle, M., 2019. Multicomponent ionic transport modeling in physically and electrostatically heterogeneous porous media with phreeqcRM coupling for geochemical reactions. Water Resour. Res. 55, 11121–11143. https://doi.org/10.1029/2019WR026373

Oko, C.O., Ndjaka, A.T., Minhas, R., Lassin, A., Azaroual, M., Crocollo, F., Bataller, H., 2024. Study of salt free-diffusion by 1d transport numerical simulations and shadow-graph experiments. J. Hydrol. 639, 131643. https://doi.org/10.1016/j.jhydrol.2024.131643

Pitzer, K.S., 1991. Activity Coefficients in Electrolyte Solutions. CRC Press.

Radko, T., 2005. What determines the structure of thermohaline staircases? J. Fluid Mech. 523, 79–98. https://doi.org/10.1017/S002211200400226X

Rard, J.A., 2021. Mutual diffusion coefficients of aqueous na2SO4 at c = 0.02989 mol/dm³ and =298.15 k by using Rayleigh interferometry with free diffusion boundary conditions: experimental test of the effect of the c³/² concentration dependence of refractive index at low concentration. J. Solution Chem. 50, 886–893. https://doi.org/10.1007/s10953-021-01092-5

Rard, J.A., Albright, J.G., Miller, D.G., Zeidler, M.E., 1996. Ternary mutual diffusion coefficients and densities of the system z1NaCl + (1-z1)Na2SO4(aq) at 298.15 k and a total molarity of 0.5000 mol dm-3. J. Chem. Soc., Faraday Trans. 92, 4187–4197.
https://doi.org/10.1039/FT9969204187

Rard, J.A., Miller, D.G., 1979. The mutual diffusion coefficients of NaCl·H₂O and CaCl₂·H₂O at 25°c from rayleigh interferometry. J. Solution Chem. 8, 701–716. https://doi.org/10.1007/BF00648776

Rehman, A.N., Bavoh, C.B., Lal, B., Sabil, K.M., Sangwai, J.S., 2022. Insights on CO2 hydrate formation and dissociation kinetics of amino acids in a brine solution. Ind. Eng. Chem. Res. 61 (37), 13863–13876. https://doi.org/10.1021/acs.iecr.2c02178

Shevtsova, V., Sechenyh, V., Nepomnyashchy, A., Legros, J.C., 2011. Analysis of the application of optical two-wavelength techniques to measurement of the soret coefficients in ternary mixtures. Philos. Mag. 91, 3498–3518. https://doi.org/10.1080/14786435.2011.586376

Steefel, C.I., Lichtner, P.C., 1994. Diffusion and reaction in rock matrix bordering a hyperalkaline fluid-filled fracture. Geochim. Cosmochim. Acta 58, 3595–3612. https://doi.org/https://doi.org/10.1016/0016-7037(94)90152-X

Tabrizinejadas, S., Carrayrou, J., Saaltink, M.W., Baalousha, H.M., Fahs, M., 2021. On the validity of the null current assumption for modeling sorptive reactive transport and electro-diffusion in porous media. Water 13 (16). https://doi.org/10.3390/w13162221

Trevelyan, P.M.J., Almarcha, C., D. E. Wit, A., 2011. Buoyancy-driven instabilities of miscible two-layer stratifications in porous media and hele-shaw cells. J. Fluid Mech. 670, 38–65.

Turner, J.S., Stommel, H., 1964. A new case of convection in the presence of combined vertical salinity and temperature gradients. Proc. Nat. Acad. Sci. 52 (1), 49–53. https://doi.org/10.1073/pnas.52.1.49

Šeta, B., 2020. Determination of Diffusion and Thermodiffusion Coefficients in Multicomponent Mixtures. Numerical and Experimental Analysis of Gravitational Convective Instabilities. Ph.D. thesis. Universitat Rovira i Virgili. Tarragona.

Šeta, B., Errarte, A., Dubert, D., J. Gavaldà, J., Bou-Ali, M.M., Ruiz, X., 2019. Gravitational stability analysis on double diffusion convection in ternary mixtures. Acta Astronaut. 160, 442–450.
https://doi.org/https://doi.org/10.1016/j.actaastro.2019.04.020

Šeta, B., Spangenberg, J., Bou-Ali, M.M., Shevtsova, V., 2025. Impact of cross and main diffusion coefficients on symmetry breaking in nonreactive diffusion systems. Phys. Rev. E 112, L033101. https://link.aps.org/doi/10.1103/yr2z-b9px. https://doi.org/10.1103/yr2z-b9px

Yadav, D., Kim, M., 2015. The onset of transient soret-driven buoyancy convection in nanoparticle suspensions with particle-concentration-dependent viscosity in a porous medium. J. Porous Media 18 (4), 369–378. https://doi.org/10.1615/JPorMedia.v18.i4.10

Yan, L., Niftaliyev, R., Voskov, D., Farajzadeh, R., 2025. Dynamics of salt precipitation at pore scale during CO₂ subsurface storage in saline aquifer. J. Colloid Interface Sci. 678, 419–430.