Experimental and Simulation Study on Validating a Numerical Model for CO2 Density-Driven Dissolution in Water

Carbon dioxide density-driven dissolution in a water-filled laboratory flume of the dimensions 60 cm length, 40 cm height, 1 cm thickness, was visualized using a pH-sensitive color indicator. We focus on atmospheric pressure conditions, like in caves where CO2 concentrations are typically higher. Varying concentrations of carbon dioxide were applied as boundary conditions at the top of the experimental setup, leading to the onset of convective fingering at differing times. The data were used to validate a numerical model implemented in the numerical simulator DuMux. The model solves the Navier–Stokes equations for density-induced water flow with concentration-dependent fluid density and a transport equation, including advective and diffusive processes for the carbon dioxide dissolved in water. The model was run in 2D, 3D, and pseudo-3D on two different grids. Without any calibration or fitting of parameters, the results of the comparison between experiment and simulation show satisfactory agreement with respect to the onset time of convective fingering, and the number and the dynamics of the fingers. Grid refinement matters, in particular, in the uppermost part where fingers develop. The 2D simulations consistently overestimated the fingering dynamics. This successful validation of the model is the prerequisite for employing it in situations with background flow and for a future study of karstification mechanisms related to CO2-induced fingering in caves.


Introduction
The dissolution of CO 2 in water slightly increases the density of water; e.g., [1]. This can lead to an unstable layering of CO 2 -rich water with a higher density above water with lower density. Convective mixing characterized by protruding fingers is triggered after a certain onset time.
In geological systems, we often see double-diffusive convective systems, where density differences typically arise from changes in the concentration of the dissolved components and from differences in temperature. For this study at atmospheric pressure, we focus on concentration-dependent processes by keeping the temperature constant, and consequently denote the process as density-driven or density-induced dissolution.
Convective dissolution is of practical importance in various technical and geological fields of research. One such field is CO 2 geological storage. The injection of supercritical CO 2 into a geological formation, e.g., a saline aquifer, typically leads to a buoyancy-driven segregation of the CO 2 and the brine. The CO 2 fluid density depends on the pressure and temperature in a reservoir and is often in the order of half the density of the brine. Thus, the CO 2 phase will end up in a stratum underneath a hydraulic barrier on top of the brine. Over time, CO 2 dissolves in the brine, increases the brine's density, and triggers a fingering process, which eventually results in an enhanced dissolution and an effective vertical downward transport of CO 2 . This effect has already been discussed in early publications in the field of CO 2 geological storage, e.g., [2,3], and is denoted also as solubility trapping [4]. On the subject of CO 2 geological storage, many publications are found on (in-)stability analyses and estimates for the time until the onset of fingering or the wave length of the fingering pattern in porous media, e.g., [5][6][7][8][9], or the scaling with different dimensionless numbers, e.g., [10]. High-resolution numerical studies on Darcy-type models for porous media also show that the spatial discretization length has to be very small relative to the scale of a typical storage reservoir in order for us to resolve onset time and fingering pattern correctly [7,11], making grid-converged results on large spatial scales practically infeasible. Therefore, more pragmatic approaches avoid the resolution of the fingers and employ effective rates, dependent on permeability; density difference as a function of CO 2 concentration and brine salinity; and fluid viscosity [11]. An overview is given, for example, in the paper of Green and Ennis-King (2018) [12].
This study is motivated by another problem, featuring similar processes of CO 2 density-driven dissolution but in a completely different context. Karstification is a field where CO 2 density-driven dissolution in a water body may play a role underestimated so far. Scherzer et al. (2017) [13] formulated the concept of nerochytic (from Greek nerochytis = sink) speleogenesis which hypothesizes the density-driven dissolution of CO 2 from cave air with elevated CO 2 concentrations in the cave atmosphere in the water body to be a relevant mechanism for karstification. CO 2 density-driven dissolution might be a process that can be considered an essential supplement to common karstification models, since it could enrich water bodies with additional carbonic acid from cave air with its seasonal fluctuations in CO 2 concentrations. Thus, it would contribute to a perpetuation of the dissolution potential in the phreatic zone and consequently increase the dissolution rate of CaCO 3 there. In this context, we are interested in estimating the time-scales on which an enrichment of resting cave waters with dissolved CO 2 may take place. Furthermore, a quantitative understanding of the seasonal dynamics of CO 2 input to cave water may allow estimates on how significant this transport process might be for karstification. We consider numerical modeling as an appropriate tool for studying such effects over the relevant geological time-scales, including both the dynamics in flow and transport processes, and the chemistry of the CO 2 reacting with water to form carbonic acid and further upon dissolving calcium carbonate. Our vision is to develop and validate such a model; this study aims at contributing to it. The particular focus of the validation attempt is on the sinking velocity of the density-induced fingers, since, for the application we have in mind, it is far more important to capture the dynamics of the fingers during seasonal fluctuations in CO 2 concentrations rather than the onset time of the first fingers. According to [14], validation is the substantiation that a computerized model within its domain of applicability possesses a satisfactory range of accuracy consistent with the intended application of the model.
In other words, we consider validation to be the successful comparison of the simulation result with data from a well-controlled experiment. In this study, we present a first step towards this goal.
We set up a quasi-two-dimensional experiment in a laboratory flume with dimensions 60 cm × 40 cm × 1 cm, where we visualized density-driven dissolution of CO 2 using a pH-sensitive color indicator; see Section 2.1. The data of onset time and evolving fingering pattern over time were then compared with the results of a numerical model that solves the Navier-Stokes equation with a density dependent on CO 2 concentration, and a transport equation for CO 2 including advective and diffusive terms. The model is presented in Section 2.2 and the results and discussion of the comparison between the experiment and simulation are found in Sections 3 and 4.

Experimental Setup and Procedure
The experimental setup is shown in Figure 1. The flume has the dimensions 60 cm × 40 cm. The metal frame structure and sealing allow an effective length of 57 cm. The effective water-filled height is 32.5 cm. The aperture between the glass plates at the front and the back is effectively 1 cm. Thus, the flow induced in the flume is quasi-2D, while the aperture is still large enough to allow for the development of a velocity profile in this direction. This is discussed later on in the modeling section, where the differences between 2D and 3D simulations are analyzed. The flume can be filled and emptied with flexible tubes to the left and right; see Figure 1a. For this study it was filled with deionized water and bromocresol green as a pH indicator, as in the study of Thomas et al. (2015) [15], who showed that this technique enables the visualization of fingering in density-driven dissolution. Bromocresol green is deep blue at pH values higher than 5.4, and changes to green, then yellow between 5.4 and 3.8. It was added to the water at a concentration of 7.5 × 10 −5 mol/L to allow for good color contrasts while not affecting the fingering dynamics too much. The color indicator consumes and buffers part of the acidity that is in the system when CO 2 dissolves and forms carbonic acid, and thus acts as a buffer and slightly slows down the fingering dynamics.  At the top of the flume, 1 cm above the water table, a very fine metal grid separated the water from a gas flow of ≈3 l/h that was induced from right to left using a squeeze pump in a flexible tube. This means the gas was flowing at a velocity of ≈0.2 cm/s, i.e., at a low Reynolds number (see Equation (14)) and causing very low pressure losses between inlet and outlet. The gas flow established a defined concentration of CO 2 as a boundary condition for the water in the top part of the flume. However, in accordance with the small flow velocity, this occurred gradually within ≈5 min from right to left. The inserted metal grid aimed at minimizing shear stresses at the water surface and perturbations in the water. The gas was fed at a fixed CO 2 concentration from a 100 L bag into the flume and in a closed circuit back into the bag. The large volume of the bag aimed at keeping the CO 2 concentration in the gas, and thus in the water table, virtually constant.
Experiments were carried out under atmospheric pressure conditions, but for different CO 2 concentrations in the gas flow, at 8 • C in a cooling chamber and a few others also at 20 • C. The temperature affects several quantities in this setup: equilibrium CO 2 concentration in water is smaller at higher temperature, water density and water viscosity both decrease with higher temperatures.
A Nikon D5000 SLR camera was used to take 2 photos per minute with aperture f/8, exposure time 1 100 s, ISO 500. The photos were processed to improve image contrast and to crop out the glass front of the flume. The green light beam was then extracted to obtain an image as in Figure 2a where each pixel is defined as an integer value between 0 (white) and 255 (black). A threshold value in the order of 120 allowed distinction between those pixels that showed a state where the color indicator had already changed, i.e., pH ≤ 5.4, and the other pixels where the change was not yet visible. In this way, a (0,1) matrix was obtained as illustrated in purple (0) and yellow (1) in Figure 2b. This illustration then served for the evaluation of the fingering pattern. For the calculation of the finger's downward-oriented velocities, the last 10% at the left and right edge of the flume images was cut off in order to put the focus only on the fingers without influence from the boundaries.

Governing Equations
The mathematical description of density-induced convection in our case involves the continuity equations for the components water and CO 2 , both present in the aqueous phase, as well as the Navier-Stokes equations. The dissolution process at the interface between the liquid and the ambient air is not modeled explicitly. Instead, the gaseous concentration of CO 2 in the ambient air is converted via Henry's Law into an aqueous CO 2 concentration in equilibrium and is accordingly assigned as a fixed value at the top boundary.
The continuity equation for each component κ ∈ {H 2 O,CO 2 } is given by is the density of the aqueous phase and dependent on the concentration of CO 2 , expressed here by the mass fraction X. The velocity vector is denoted by v and D is the binary diffusion coefficient.
The Navier-Stokes equation is formulated here as µ is the dynamic viscosity and depends on temperature, p is pressure, and g is the gravitational acceleration vector.
For the sake of computational efficiency, flat three-dimensional domains may be simulated with a two-dimensional numerical model. Assuming a parabolic velocity profile along the axis of the omitted dimension, Flekkøy et al. (1995) [16] proposed a friction term which is added to the right-hand side of Equation (2). This "pseudo-3D" approach has been used successfully for different applications, provided that the depth in the neglected dimension is sufficiently small compared with the other dimensions of the domain [17][18][19]. The constant c determines whether the maximum velocity in the central plane at 0.5h (c = 8) or the height-averaged velocity is considered (c = 12). In this work, we provide numerical results calculated in full 3D, 2D, and pseudo-3D using c = 12.

Numerical Solution
The numerical simulator DuMu x (www.dumux.org) provides the platform for solving the system of equations. All implementations used for this study can be reproduced and found for download at https://git.iws.uni-stuttgart.de/dumux-pub/class2019a.
Here, we use the freeflow Navier-Stokes model in DuMu x and employ the brineco2 fluid system. Two scalar values, i.e., pressure and concentration/mole fraction, and the velocity vector, are selected as primary unknowns for which the system of equations is solved. We use a staggered-grid method that corresponds to a finite-volume method with different control volumes for different equations. The control volumes for the velocity components are shifted from the control volumes for the pressure and mass fractions. The staggered-grid method is a robust scheme without pressure oscillations and mass conservation. All equations are solved fully implicitly in time, using a Newton method to treat non-linearities. The Newton scheme chooses the time-step size according to its convergence with a user-controlled maximum time-step size. For further details on discretization, numerical solution methods, and their implementation, we refer to [20] or the handbook of DuMu x [21].
We provide a brief summary of basic model assumptions as far as they are important for this study. Justified by slow processes of diffusion and density-driven advection, we assume equilibrium between the gaseous and aqueous CO 2 concentration at the gas-water interface. This allows us to assign the aqueous CO 2 concentration as a boundary condition at the top. We compute both the density and viscosity of water as functions of temperature and pressure, while the density also depends on the dissolved CO 2 concentration. The binary diffusion coefficient for CO 2 in water is assumed to be constant, since its variation with temperature, pressure, and concentration is very small [22,23].
Three different setups for the numerical simulations were employed; namely, a 2D model, a 3D model, and a 2D model with a correction term for wall friction which we denote "pseudo-3D." Dirichlet conditions were assigned to the top boundary, where the pressure was fixed at atmospheric conditions, the velocity components were set to zero, and the mole fraction of CO 2 was calculated via Henry's Law from the given partial pressure of CO 2 ; see Equation (7). All other boundaries were no-flow/no-slip boundaries. The numerical grids used for this study are explained further below. The maximum time-step size in the simulations was limited to ∆t max = 30 s.

Grids Used in The Simulations
Previous studies in the literature, e.g., [7,11], showed that numerical grids need to be highly refined in order to capture the onset time and the dynamics of the fingers. In this study, we also tested the influences of grid resolution on these phenomena and compared the results with the experimental observations. We started with a regular grid and a discretization length of δx = δz = 0.006 m both in 2D and in 3D, and δy = 0.0017 m in the third dimension of the 3D model. Accordingly, the grids had 5280 nodes in 2D and 36960 in 3D respectively. For comparison, we refined the grid in the topmost part of the flume and applied a grading which resulted in very small elements at the top. Figure 3 shows the graded grid in full (left figure) and in a zoom into some of the topmost cells (right figure) where the grading is applied. Underneath the graded region, both grids, the regular and the graded, have the same spatial discretization of 0.006 m. The smallest grid cell at the top of the graded region has 0.0003 m in the vertical extent. The grading is applied since the region at the top of the water body, i.e., close to the gas-water interface, is where diffusion takes place and where the essential parameters that characterize instability (see below) need high accuracy. This refers mainly to a high spatial resolution of CO 2 concentrations due to diffusion, and, corresponding to the CO 2 concentration, the density variation which eventually drives the fingering dynamics.

Density Variation
The partial differential equations are coupled via the density which depends on the CO 2 concentration. We follow here an approach suggested by Garcia (2001) [1].
The density (in kg/m 3 ) is computed as H 2 O is the density of pure water dependent on pressure and temperature, M H 2 O is the molar mass (in kg/mol) of pure water, and M T is accordingly obtained from The apparent molar volume of dissolved CO 2 V φ (in m 3 /mol) is calculated as a function of temperature T (in • C) from Reference pure-water densities in this study were 999.85 kg/m 3 at 8 • C and 998.21 kg/m 3 at 20 • C, as implemented in our simulator based on IAPWS formulations [24].

Henry's Law for Calculating Aqueous CO 2 Concentration
The concentration of CO 2 at the interface between the atmosphere and the water body is calculated as a function of the partial pressure of CO 2 p CO 2 (in atm) in the ambient atmosphere by assuming equilibrium between the fluid phases. Accordingly, Henry's Law is assumed to be valid: H aq,CO 2 (in mol CO 2 /mol H 2 O·atm) is the temperature-dependent Henry constant for CO 2 in water.
For this study, we evaluate Henry's Law at two different temperatures in our experiments. Using reference values and corresponding equations of state from the literature [25], we obtain values for Henry's constant as 9.99 × 10 −4 mol CO 2 /mol H 2 O·atm at 8 • C and 7.04 × 10 −4 mol CO 2 /mol H 2 O·atm at 20 • C.

Characterizing Instability
Density-driven convection is a complex phenomenon that occurs when a fluid of higher density rests upon a fluid of lower density. It was described phenomenologically and mathematically by Bénard (1901) [26] and Rayleigh (1916) [27] for convective cells forming due to density differences induced by a fluid of lower temperature resting upon a fluid of higher temperature. Our focus is on the enhanced dissolution of CO 2 upon the onset of finger-like swathes of water with more dissolved CO 2 moving downward, and others with lower concentration moving upward while fostering more dissolution at the liquid-gas interface. Density differences are caused solely by varying CO 2 concentrations while the temperature is kept constant at 8 • C or 20 • C.
The dimensionless Rayleigh number is commonly employed to characterize instability. It can be interpreted as the ratio of a characteristic diffusion time to a characteristic convection time. Important factors of influence are the density difference ∆ , the diffusion coefficient D, the fluid's dynamic viscosity µ, and the spatial dimensions of the setup. The setup of our study has a thickness d of only 1 cm, resembling a Hele-Shaw type of cell. We therefore use a definition of the Rayleigh number as it is used in porous-media applications [12] with the permeability k (in m 2 ) representing a resistance to advection and a porosity equal to one. For flow between parallel plates, the permeability can be derived as k = b 2 /12, where b is the distance between the plates; in our case, 1 cm. Using H as the vertical dimension of the setup, i.e., 40 cm, Rayleigh's number is then formulated as High Rayleigh numbers are favorable for convection, while diffusion is dominant in low Ra regimes. The Rayleigh number is thus a criterion for the onset time of a convective fingering and also for the characteristic wave length of the fingers. Again following Green and Ennis-King (2018) [12], Pau et al. (2010) [11] or Emami-Meybodi et al. (2015) [28], who review the literature on convective dissolution of CO 2 in saline aquifers, we may estimate for the onset time based on linear stability analysis: c 0 is a proportionality factor. Similarly, the critical wave length λ crit can be approximated by where c 1 is again a proportionality factor. As we want to validate a numerical model for a physical process where density-driven convection competes with diffusion, we consider the dimensionless Péclet number an important criterion. Here, we use a Péclet number with the spatial discretization length ∆l as the characteristic length and formulate v c is the characteristic convective flow velocity, which consequently depends on the concentration-induced density difference. Since we use a porous-media analog here, we estimate v c with a Darcy-type approach as It is well known that Pe < 2 is a theoretical criterion for numerical stability in a linear convection-diffusion problem; e.g., [29]. Ensuring Pe < 2 guarantees that the dynamics of the processes can be captured properly on a given grid. We can argue similarly with respect to the temporal discretization. We note that the numerical scheme employed in this study (upwinding for flux terms, fully implicit in time) does not require criteria for numerical stability. The Courant number which is restricted to Cr < 1 for explicit time discretization, allows us to define a maximum time-step size ∆t for an optimal resolution of the physical processes in time for a given grid-cell size. Table 1 gives an overview of the metrics we introduced for characterizing finger instability. Note that we have performed simulations and experiments at two different temperatures and different CO 2 partial pressures in the gas above the water table. The values indicate that higher instabilities should be expected at lower temperatures. The "optimal" values for discretization length and time-step size are obviously very small and totally unfeasible for numerical simulations. We applied 6 mm as the spatial discretization length in the regular grid and in the graded grid below the graded region. The maximum time-step size was restricted to 30 s. A grid convergence study with discretization lengths of 1.5 cm, 1 cm, and 6 mm in a regular grid showed that the pattern of the fingers did not change significantly in the step from 1 cm to 6 mm, while onset times were still varying; this motivated our decision to use the graded grid. Given the uncertainty associated with the buffering capacity of the color indicator, we hold the numerical accuracy achieved with this choice as reasonable.
Let us finally define also the Reynolds number as

Figures 4 and 5
show the comparison between pseudo-3D numerical simulation results, obtained with the graded grid, and pictures taken at three different times during the experiments with p CO 2 = 0.25 atm ( Figure 4) and with p CO 2 = 0.75 atm ( Figure 5), both at 8 • C. The pseudo-3D approach is computationally as efficient as the 2D approach, and in terms of the relevant metrics, i.e., onset time and finger velocities, almost as "correct" as the 3D approach, as shown further below. Numerical results are displayed with a constant color legend for the mole fractions in all plots within each figure, but minimum/maximum values of CO 2 mole fraction differ from figure to figure. We did not attempt to derive reliable values of CO 2 concentrations from the photos, since this would have required a very costly calibration procedure, probably without much additional benefit for this study. The color indicator consumes some CO 2 before changing color. Therefore, we do not consider it appropriate to try to compare the very details between the numerical results and the photos. The color indicator might suppress the diffusive part of the processes, whereas we expect that the dynamics of the convective fingers are reasonable to compare. This is indeed confirmed by the comparison provided in Figure 4 and 5 with a very satisfactory agreement between model results and experiments.
Experiments at 8 • C were successfully performed for 1 atm, 0.75 atm, 0.5 atm, 0.25 atm, and 0.05 atm. We also conducted two experiments at 20 • C at p CO 2 = 1 atm, 0.75 atm, and 0.5 atm. Without any fitting of parameters, the onset of instability is at about the same time in model and experiment when the graded grid is used. The same holds for the progress of the fingers. At later times, there is a tendency towards more disagreement between numerical results and experiments. The experimental fingers seem to get stuck or retarded after some time, and cease to show the same dynamics as the numerical fingers. The buffering by the color indicator probably hinders the visualization of the diffusive processes that tend to dilute the protruding fingers as they grow.
A full set of photos is provided in the Supplementary Materials; it is representative of the experimental setups with 0.75 atm both at 8 • C and at 20 • C.
The initial dynamics are mainly determined by the evolving convective flow once the instability sets in. Therefore, as a criterion for model validation, we will compare the simulations and experiments with regard to the onset times and the velocities of the protruding fingers.  The 2D simulations neglect the third dimension, while the pseudo-3D model accounts for it phenomenologically by a drag term for wall friction. In fact, the 1 cm distance between that back and front plate has an influence on the dynamics, since no-slip conditions slow down the induced flow. Using k = b 2 /12 to estimate a permeability analogous to that in a porous porous medium yields a value of 8.3 × 10 −6 m 2 . This is a relatively high value, and a Darcy model is therefore not appropriate since Reynolds numbers at typical finger velocities in this study are in the order of 1 or slightly higher; i.e., too high. Figure 6 shows an exemplary snap-shot of a 3D simulation with p CO 2 = 1 atm at 8 • C. It is clearly visible that the CO 2 concentrations accumulate at all no-flow boundaries. Thus, density differences are high there in particular, and the onset of fingering starts at the boundaries. This effect is not seen in the 2D and pseudo-3D simulations, while the images of the experiments confirm the accumulation of CO 2 at the lateral boundaries. The grid resolution in the simulations is very high in the top few centimeters, and thus almost comparable with high-resolution studies; for example, [7,11]. We tested further grid refinement below the graded region only for the 2D simulations, but this did not further improve the mainly qualitative comparison with the experiments. Onset times were determined in both simulation and experiment simply as those times when the initiation of fingers could be clearly discerned in the plots and images. Of course, this involves some uncertainty, but the overall significant uncertainties in the experiment justify this approach. Similarly, the wave lengths were "measured" from the finger tips observed in the simulations. We did not try to express the wave lengths in numerical values. A qualitative look at the images and comparison with the simulations shows similar behavior; see also Figures 4 and 5. The onset times of the pseudo-3D simulations are summarized in Table 2. 3D results and their evaluation were very similar and are omitted for the sake of brevity. The table further provides the values for the proportionality factor c 0 as obtained by solving Equation (9). The values are not constant, but rather differ over orders of magnitude; thus, they clearly indicate that the factors of influence in Equation (9) do not cover all relevant processes properly.   Figure 7 shows the onset times depending on the concentration of CO 2 in the atmosphere at the top boundary and depending on the Rayleigh number. The plots in Figure 7 use model results from the graded grid. They show a clear and consistent tendency towards later onset for smaller CO 2 concentrations and smaller Rayleigh numbers. We can see the same tendency in the plots of Figure 8, which uses the model results from the regular grid. The experimental onset times are in very good agreement with the simulations on the graded grid, except for the case with 0.05 atm where the experimental value at 8 • C is earlier than the models'. At 8 • C and 0.5 atm, for example, we see deviations from the experimental value of 0.8% higher for the 2D graded grid, 18% lower for the 3D graded, 7.5% higher for the pseudo-3D graded, 88% higher for the 3D regular, 33% higher for the 2D regular, and 99% higher for the pseudo-3D regular. Experimental results at 20 • C are not available for all CO 2 concentrations. While the results, both in the simulations and the experiments, did not show a strong influence of the temperature, we faced problems visualizing the fingers with the color indicator for higher temperatures and smaller CO 2 concentrations. The results on the regular grid do not match as well as those on the graded grid. Here, the 2D results seem to match better than the pseudo-3D and 3D. We argued above that pseudo-3D and 3D should be expected to be more realistic, since they consider wall friction. Therefore, further explanation is required to understand this effect. In fact, this depends on the implementation of the Dirichlet boundary condition for the CO 2 concentration at the top. Initially, we assumed that no dissolved CO 2 was present in the topmost cell. Therefore, at time zero, diffusion starts filling up the first cell and instability develops in terms of a growing density difference. The smaller the topmost cell, the faster the concentration in this cell can increase, which essentially leads in a coarser grid to an underestimated concentration, and thus, underestimated density. Generally, a coarse grid then shows a tendency towards overestimating the onset time, which is partly "compensated" by the neglected wall friction in the 2D model. This explains the seemingly better matching 2D model in the regular grid.  We consider the velocities of the protruding fingers an important criterion in the comparison between experiments and simulations, since they describe the dynamics of such a setup in the long run. The finger velocities were determined from the images, as explained in Section 2.1; see also Figure 2. The same approach as with the images from the experiment has been applied to the plots from the simulations. We used a range between 1 × 10 −6 and 3 × 10 −6 for the mole fraction of CO 2 in water to identify fingers with sufficiently sharp contrasts. We provide the details in Supplementary Materials to this manuscript. With this procedure, the characteristic velocities were determined; they are plotted in Figure 9 for the graded grid and in Figure 10 for the regular grid, in both cases depending on the CO 2 partial pressure and on the Rayleigh number respectively. All plots reveal a tendency towards faster fingers at higher concentrations, and thus higher density differences. Using the same example as for the onset times, i.e., at 8 • C and 0.5 atm, we observed deviations from the experimental value of 393% higher for the 2D graded grid, 25% higher for the 3D graded, 73% higher for the pseudo-3D graded, 33% lower for the 3D regular, 190% higher for the 2D regular, and 11% lower for the pseudo-3D regular. The first obvious observation is that the finger velocities from the 2D model on the graded grid are far from all other data points. This is proof of the importance of considering the wall friction, as we do in the pseudo-3D and in the 3D model. On the regular grid, it is also the 2D finger velocities which have the highest values, but they are much closer to all other data points than on the graded grid. This, in turn, confirms our argument from above that the density difference in the top layers is underestimated on the coarser grid. For the pseudo-3D and the 3D cases, the grid does not have such an obvious effect in the comparison with the experimental data, but with the explanations above, we have clearly made our case for the pseudo-3D model as the best trade-off between desired physical accuracy and computational efficiency.  The "measurement" of fingering velocities is prone to be affected by random effects in the fingering patterns. This might to some extent explain the noise in the data and some outliers we see, for example, at 0.75 and 1 atm, which, however, seem less pronounced in the plot over the Rayleigh number.
Wave lengths were evaluated only qualitatively. They were slightly larger for smaller CO 2 concentrations, as expected, but were not analyzed further.

Discussion
The results show a very good agreement between experimental results and numerical simulations on the graded grid, and they were still reasonable on the regular grid. The metrics we used for the comparisons were rather simple and qualitative, but we consider this appropriate given the uncertainties in the experiment. The color indicator has proven useful to visualize the dynamics of the convective fingers, while it compromised the diffusive effects that would counteract the convective growth of the fingers by diluting the concentrations at the finger tips. The qualitative comparison of the density-driven concentration pattern as shown in Figures 4 and 5 looks promising, although for later times, the match between experiment and simulation is impaired mainly by the above-mentioned limitations in the experiment. In particular, we note that the shift in onset time and temporal evolution from Figure 4 to Figure 5 is in good agreement between experiment and model. The qualitative and quantitative difference in the finger's dynamics from one scenario to the next is in very good agreement regarding experimental and numerical results. Thus, this study is a significant step towards a successful validation of our Navier-Stokes model with concentration-dependent density.
The 2D simulations neglect the no-slip conditions at the front and back plate of the flume. Thus, as expected, they tend to overestimate the finger growth and propagation velocity, but, as the comparison with the experiment and the 3D simulations shows, this is still within the same order of magnitude. Both the 3D and the pseudo-3D results show a close match in terms of reproducing the characteristic velocity of the fastest fingers; see Figure 10. Our setup with 1 cm distance between front and back plate is a trade-off between our goals (i) to have the processes in an open water body, and (ii) to use simple imaging by taking a photo of the front glass plate as representative of the fingers developing inside.
Onset times are consistently overestimated in the simulations on the regular grid compared with the experiment. They match well when the graded grid is used, except for the case of a very small CO 2 concentration being present, where the experimental onset is seen significantly earlier than in all models. The experiment is probably disturbed by non-perfect conditions at the top boundary. The constant gas concentration at the top was achieved not instantaneously, but gradually, according to the velocity of a gas flow along the top of the water body, which was, however, separated from the water by a metal grid to minimize disturbance. Nevertheless, even if small, such a disturbance might lead to a premature triggering of fingers. If so, it does not seem to be very important for a large enough CO 2 concentration, but for small concentrations where onset is late, it seems to play a role. On the other hand, the buffering by the color indicator should act protractingly. A third point: as we have shown, the implementation of the discretization scheme (see Section 2.2.2) [20,21] can cause a delay of the onset time for coarse grids, since it requires some time to have the first discrete finite volume at the top of the water body at a high enough CO 2 concentration to trigger instability. While it is not easy to separate these effects quantitatively, we can still state that the agreement between experiment and simulations is very reasonable and the tendency towards longer onset times for lower CO 2 concentrations is consistently reproduced in all cases.
The influence of temperature is not very distinct and clear. Since the Rayleigh numbers are dependent on the temperature, see Table 1, the instability at the same CO 2 concentration should be higher at 8 • C than at 20 • C. Consequently, this would favor an earlier onset and faster fingers at the lower temperature. On the one hand, the density difference is influenced by temperature due to the solubility being temperature-dependent. On the other hand, there is the viscosity as the dampening effect for developing fingers, which shows the opposite tendency. Viscosity is high at low temperatures; thus, the higher instability encounters a more viscous fluid and vice versa. While at 8 • C, the viscosity is higher by a factor of 1.35 than at 20 • C, the density difference at 8 • C is higher by a factor of 1.41 than at 20 • C. Both parameters occur linearly in the numerator and denominator of the Rayleigh number. Therefore, since the relative change is very similar, it is expected that in this case, this has no significant impact on the onset times. Regarding the finger's velocities, it is not as simple, since the density difference diminishes during the advancement of a finger and viscosity might become the stronger argument. In any case, as we expected, the data from the experiments and models do not give us strong evidence for a temperature dependence.
With respect to the choice of the model, we can state that the 2D model consistently overestimates the dynamics of the fingers, and does so very strongly on the graded grid. While the impact on the onset time is visible, it is much more important for the velocity of the fingers. Shear stress is imposed by the back and front glass plate when fingers start developing. The pseudo-3D model accounts for this with a correction term in the momentum balance. Regarding the metrics we apply, i.e., onset time and velocity of the fingers, the pseudo-3D model achieves very good agreement, since it is close to the results of the full 3D model in all instances. We consider the velocities of the fingers as the most important metric for comparison, and here we clearly see that pseudo-3D and 3D both work well, while the neglection of wall friction in the 2D model is the main reason for the deviations observed.
With regard to the characterization of instability, we tried to use metrics from the abundant literature on the convective dissolution of CO 2 injected into saline aquifers, where Equations (9) and (10) can be used to predict onset time and critical wave length. In our study, we assessed the value of c 0 for the onset time, and found that this proportionality factor is far from being a constant value. Interestingly, when c 0 is divided by the Péclet number as given by Equation (11), the values obtained are almost constant. This is shown in Table 3  While the values for c 0 vary by a factor of about 60 between the respective minimum and maximum values (Table 2), this variation is reduced to a factor of about 3 for c 0 /Pe (Table 3). Using the definition of the Pe number from Equation (11), we can write From Equation (15), one can conclude that, with g, k, and ∆l all being constant, and µ varying only slightly with temperature, an almost constant c 0 /Pe would suggest a rather anti-proportional correlation between t onset and ∆ . In contrast to that, Equation (9) suggests that ∆ has a quadratic influence on t onset .
We note that we did not pursue further theoretical considerations or stability analyses on the given system of equations consisting of Navier-Stokes equation and mass balance equations for the two components coupled via a density depending on concentration. Our study is rather descriptive based on observations, which we have found to be in satisfactory agreement between numerical simulations and experiments.
Finally, let us offer some remarks on the computational efforts in the three different models. We may compare the CPU time required to simulate 2000 s in the case of 8 • C with 50 % CO 2 as the boundary condition using the graded grid. The full 3D model required 152,800 s, the 2D model 525 s, and the pseudo-3D model 325 s. Obviously, the 3D model is by far the most costly. In fact, the pseudo-3D shows the best performance, even in comparison to the 2D model. The tiny additional effort of computing the drag term is easily compensated for in this case by a better convergence of the Newton solver, since velocities in the pseudo-3D are smaller than in the 2D simulations.

•
This study is a significant step towards a successful validation of the applied numerical model. It shows very good agreement for the onset times, for the characteristic finger velocities, and for the qualitative fingering patterns between experiments and simulations.

•
We aim at applying this model concept in future studies for the evaluation of the dynamics of CO 2 migration in the phreatic zone of caves during seasonal fluctuations of CO 2 concentrations in cave air. In order to evaluate the relevance of nerochytic speleogenesis for karstification, various aspects need to be extended. This includes geochemical processes at the water-karst interface, very complex scenarios with varying boundary conditions with seasonal variations of CO 2 concentration, and water background flow. It is therefore important to conclude from this study that the basic hydraulic processes and dynamics are well reproduced.

•
The pseudo-3D model is computationally as cheap as the 2D model, while being capable of capturing the important effect of shear stress such that the metrics for comparison show a very good agreement with the full 3D simulations. • The 2D model without the consideration of wall friction is not appropriate for capturing the fingering dynamics.

•
Our experiments provided very nice results for lower temperatures, but at room temperature we had difficulties achieving good quality images. With the concentration of the color indicator as applied in this study, the experiments at room temperature were not successful for low CO 2 concentrations.

•
Results from the literature on linear stability analyses for convective dissolution of CO 2 in saline aquifers cannot be transferred exactly to this setup. Here, we have Reynolds numbers in the order of 1 or slightly more, meaning that the requirements for a Darcy regime are violated, although not by much. In this transition regime, we observed that the relevant parameters for the onset time of fingers enter a functional description in a different order than in the relation found throughout the porous-media literature. A detailed analysis is, however, beyond the scope of this study.