# 6.3 Predictability of Convective Flow Patterns

As the example of the Elder problem showed, numerical solutions can be plagued by bifurcations. Bifurcation means that the convection regime can transition to a different steady state. Transient oscillatory behavior means no steady state exists; the flow and transport field keep changing and thus go in and out of different modes. The flow and transport field does not settle into one steady solution. Only modest hydraulic conductivity and density differences are required to achieve a Rayleigh number in excess of what is required for transient oscillatory behavior (*Ra*_{c} ≈ 240-300). Rayleigh numbers in many practical groundwater problems of interest are much higher. Post and Kooi (2003) reported values on the order of 10^{5} < *Ra* < 10^{6} for coastal aquifers in the Netherlands and van Dam and others (2009) estimated that the Sabkha aquifer in Abu Dhabi has *Ra* > 10^{4}*.*

The complex physics associated with convective processes raises important questions about the predictive ability of numerical models:

- What, if anything, can reasonably be predicted about fingering processes associated with free convection?
- How uncertain are those predictions likely to be?
- What features are physical reality (bifurcations, oscillatory behavior) and what are numerical artifacts?

Figure 22 provides preliminary insight regarding the answers to these questions. It shows the results of a laboratory tank experiment that was originally published by Post and Simmons (2010). In the original photographs of the experiment, showing the salinity distribution at different times, the view of the left part of the tank has been replaced by the results of the numerical model to facilitate assessment of the similarities and differences between numerical model results and physical systems. An accompanying video, shows the time-lapse sequence that provides more detail. The experiment was conducted to investigate free convection in the presence of heterogeneity, in this case, a low-permeability layer embedded into more permeable sand. Initially, the tank contained fresh water. As salt water was introduced from the source box at the top, then free convection occurred.

The hypothesis prior to conducting the experiment was that once the salt water reached the top of the low-permeability layer, fingers would develop in this layer but with a larger wavelength than the fingers that formed in the sand. In the laboratory experiment, this “re-fingering” did indeed occur. Interestingly, however, in the end, the low-permeability layer became salinized from below, instead of above where the salt source was located. The reason is that fresh water became entrapped below the low-permeability layer and the fresh water’s positive buoyancy caused it to drift upward, thereby pushing out the newly formed fingers and drawing in salt water from below. The video also shows numerous small plumes of fresh water escaping through the top of the low-permeability layer, rapidly mixing with the salt water.

Features observed in the tank that the numerical model replicated included: onset of free convection in the sand layer; horizontal spreading of the salt water along the top of the low-permeability feature; further descent of the salt water and its spreading along the bottom of the tank; and, upward drift of the fresh water and the salt water entering the low-permeability feature from below. Given the challenges with modeling free convection as listed above, this is a very encouraging result. At the same time, a comparison between the model and the laboratory results highlights some important differences. For example, the intricate smaller scale fingering pattern in the model was not observable in the tank. Also, the fingers that initially formed in the low-permeability layer in the tank, did not form in the model. These results illustrate that one can expect to have some success in simulating the overall behavior of free convective flow systems, but not the intricate details. This has also become apparent from some classical studies on convective flow phenomena. Combarnous and Bories (1975) showed that despite variability and complexity in the precise details of the convection cells themselves, similarity exists in geometrical structures, bifurcations, and transient oscillatory regimes as a function of Rayleigh number and slope of a layer (Figure 23). This suggests that some variables (e.g., approximate number of fingers, average speed of descent, and total mass of a plume) lend themselves more readily to prediction whereas exact number and locations of fingers are not easily predicted. Further to this, the data compilation of Cheng (1979) as shown in Figure 17 revealed that the convective flux (as indicated by the Nusselt number) is a function of the Rayleigh number of the system regardless of the details of the structure of the convection pattern itself.

# Homogeneous Systems

Based on the foregoing discussion, questions arise regarding what one can and cannot expect to be able to simulate and predict using models. Xie and others (2012) examined the predictability and uncertainty of free convection processes in homogeneous systems by considering quantitative indicators to represent plume behavior. These included the smaller-scale details of the fingering patterns themselves (e.g., number of fingers, deepest plume front) and the larger-scale features (e.g., vertical center of solute mass, total solute mass, and solute flux through the source zone). Many stochastic realizations, with the same fluid and solid properties but with random perturbations of the source concentration, were run to examine the plume dynamics and the range of quantitative indicators.

Figure 24 shows two such realizations selected randomly from the set of statistically equivalent realizations tested. Visual inspection reveals that the precise finger locations vary in space and time. The number of fingers themselves is different, but not very different. The bottom tips of the fingers at any given time appear, on average, to have descended to a similar level. In other words, there are features that appear to be reasonably consistent between the realizations. These are not entirely random and non-reproducible results so there is some level of predictive capacity that can be afforded in the analysis free convection behavior.

This is further illustrated in Figure 25, which shows the evolution of statistical features (i.e., mean and standard deviation) of four diagnostics with time for all of the stochastic realizations. These are “microscopic indicators” (i.e., number of fingers, NOF, and depth of penetration of the deepest finger, DPF) and “macroscopic indicators” (i.e., center of mass, COM, and the Sherwood number, *Sh*). The maximum and minimum values of these diagnostic variables are also plotted for comparison. All indicators display some level of variability. This confirms previous studies (Bachmat and Elrick, 1970; Mazzia et al., 2001; Schincariol and Schwartz, 1990; and van Reeuwijk et al., 2009) that exact reproducibility is neither possible nor should be expected.

**Figure ****25** – The evolution of statistical features (mean μ and standard deviation σ) of four measurable diagnostics with time for all of the stochastic realizations. These are “microscopic indicators” number of fingers (NOF) and depth of penetration of the deepest finger (DPF) and “macroscopic indicators” center of mass (COM) and the Sherwood number (*Sh*). The maximum and minimum values of all diagnostic variables are also plotted for comparison. Figure and figure caption modified from Xie and others (2012). Reprinted from Water Resources Research, 48(2), Y. Xie and others, Prediction and uncertainty of free convection phenomena in porous media, Copyright (2012), with permission from John Wiley and Sons.

The Rayleigh number in this system is 3.4 × 10^{5}, which is well into the transient oscillatory convection regime. Nevertheless, the variability demonstrated is surprisingly small with a coefficient of variation, *CV*, less than 0.31 for all the variables considered. The *CV* is highest for NOF ranging from 0.1 to 0.3 overtime, while the *CV*‘s for DPF, COM, and *Sh*, hover around 0.1, confirming that the uncertainty associated with microscopic indicators is greater than the uncertainty associated with the macroscopic indicators.

In summary, certain plume features are predictable despite the apparent randomness associated with the semi-chaotic fingering processes, at least in homogeneous systems. One should not expect to match the precise locations of fingers in space and time but predicting the number of fingers or their depth is feasible. Even greater confidence can be placed on integrated variables such as the center of mass or the flux. Rather than taking a deterministic view, free convection processes should be treated stochastically which is consistent with their fundamental physical nature. Predictions and comparisons of statistical properties based on multiple stochastic realizations must form the fundamental basis for prediction and uncertainty analysis. This shift from deterministic thinking to stochastic thinking suggests a paradigm shift is needed in the way free convection processes in groundwater systems are conceptualized, measured, modeled, and predicted.

# Heterogeneity

Discussion in the previous section focused on homogeneous systems. Spatial heterogeneity, however, exerts a ubiquitous and critical control on groundwater flow and solute transport. Many of the classical papers on free convection were developed under highly idealized conditions and more often than not conducted in homogeneous layers. Whilst some studies examined layering in geological systems and anisotropy (e.g., Nield and Bejan, 2006) these were fairly simple in comparison to the complex geologic conditions encountered in groundwater systems. However, insights from early work suggested that geologic heterogeneity and anisotropy are important controls on the conditions for onset of instability, as well as growth and decay of density-driven flow patterns. It is therefore critical to understand how the realistic geologic heterogeneity encountered in groundwater systems may influence free convection processes. Regardless of whether the geological medium consists of mixtures of sand and clay, fractured rock, or karstic limestone, the heterogeneity of the material is of fundamental importance.

Various papers have studied the effect of heterogeneity on unstable plume development (Schincariol and Schwartz (1990); Simmons et al. (2001); Prasad and Simmons (2003); Simmons et al. (2008); Nield and Simmons (2007); Simmons et al., 2010; Niederau et al., 2019). Figure 26 summarizes some of the key types of geologic structures that have been considered. These include “traditional” homogeneous materials, trending and continuous sedimentary structures, continuous vertical conduits, discrete vertical fractures, orthogonal fracture networks, and random distributions of fracture networks (of variable aperture and spacing).

Some early insights into the problems associated with applying Rayleigh numbers to predict the onset of instability in heterogeneous systems were provided by comparing model results for statistically equivalent permeability realizations (i.e., sampled from a log-normal distribution with the same mean and standard deviation) as shown in Figure 27. The top row of Figure 27 illustrates how realizations generated from the same distribution can vary, with low permeability regions shown as dark blue and high permeability regions as light green to yellow. The next four rows show simulation results after 700 years first for a homogeneous case and then for 3 cases with increasing heterogeneity. All the systems have the same Rayleigh number *Ra* = 490. A solute boundary condition is defined at the top of the model with a 5 kg m^{−3} density difference between the boundary and the background, ambient, groundwater. Instabilities were triggered by adding noise to the source concentration along the top model boundary (concentrations in each cell were randomly perturbed by up to +/−0.01 percent). The salt finger pattern differs for the four homogeneous simulations due to the random noise added to the source, but the penetration depth and finger number are comparable.

**Figure ****27** – Top row: Four realizations of a random permeability field. Second row: finger patterns for the homogeneous case. Third to fifth row: finger patterns for standard deviation of the log *k* field of *σ*=0.1, 0.2 and 0.5. For the heterogeneous simulations, the permeability field is overlaid in transparent black and white so that less permeable parts can be discerned by slightly darker zones, and more permeable parts by lighter zones. An animated version of this image is provided here.

As the variance of the permeability field is increased, the role of heterogeneity becomes clear. The lower three rows in Figure 27 show the simulation result for a standard deviation of the log *k *field of *σ* = 0.1, 0.2 and 0.5, respectively. Differences are distinct for even variability as low as *σ* = 0.1, where salt fingers preferentially occur in more permeable zones and fingering is suppressed where low-permeability zones are located near the source. As the standard deviation increases to *σ* = 0.5, the geological heterogeneity essentially dictates the fingering pattern. Instead of many small fingers, much larger structures have formed inside the more permeable parts of the aquifer. They have also traveled much deeper than the fingers in the more homogeneous simulations because the maximum permeability is higher in the zone occupied by the salt fingers than the mean permeability, allowing the salt water to sink faster. These simulation results are consistent with the findings of Prasad and Simmons (2003), who studied a stochastic version of the Elder problem and found that increased heterogeneity tends to reduce the number of fingers that form.

Not shown in Figure 27, but visible in the animation (video), is that the onset of free convection occurs after approximately four years for *σ* = 0.5 while for the homogeneous case, salt fingers do not become perceptible until 350 years. So, on the one hand, heterogeneity promotes fingering by providing the perturbations needed for instabilities to form. On the other hand, it reduces the number of fingers by channeling the salt water into more permeable zones of the aquifer. Also, where zones of low permeability are located near the solute source, the onset of free convection is delayed or prevented altogether.

It is evident that geologic heterogeneity exerts an important control on unstable plumes. Converting the *k* values to *K*_{f} in cm s^{−1} and calculating the variance of the natural logarithm allows comparing the degree of heterogeneity of the simulations of Figure 27 to values reported in the literature for well-investigated sedimentary aquifers. The fields considered here have = 0.05, 0.21 and 1.32, respectively. The relatively homogeneous sediments that make up the well-characterized Borden aquifer have between 0.24 and 0.37 (Sudicky, 1986; Woodbury and Sudicky, 1991) while the sand and gravel aquifer at the Cape Cod test site has = 0.26 (LeBlanc et al., 1991). This is on the same order as the values for the realizations in the third row of Figure 27, in which the heterogeneity has an important influence on finger development. It follows that even in relatively homogeneous aquifers, free convection is likely to be influenced by the local-scale variability of the permeability field. The more heterogeneous sediments of the MADE site (Rehfeldt et al., 1992) have = 4.5, which is much higher than the highest value considered when creating Figure 27. Under those conditions, sedimentary heterogeneity may be expected to dominate the free convection processes entirely.

A number of studies have examined free convection in fractured media (e.g., Graf and Therrien, 2007a, 2005; Murphy, 1979; Sharp and Shi, 2009; Shikaze et al., 1998; Simmons et al., 2008). Graf and Therrien (2007b) observed that concentration distributions were very sensitive to the number of equidistantly distributed fractures near the source, and that fracture connection pathways and fracture permeabilities controlled the depth to which the dense plume ultimately migrated. Simmons and others (2008) studied the possible modes of free convection in fractured media. Inter-fracture convection (Mode 1), intra-fracture convection perpendicular to the face of the fracture (Mode 2A) and intra-fracture convection parallel to the face of the fracture (Mode 2B) were considered. These are illustrated in Figure 28.

By developing and applying Rayleigh numbers for each possible mode of convection, the likelihood of each mode was determined. All modes of convection are theoretically possible in fractured rock aquifers at extremely small density contrasts. Mode 2B is the most likely mode followed by Mode 1 and then Mode 2A. Three-dimensional simulation is required to permit all convection modes to occur but is computationally expensive. The findings have consequences for the simulation of convection in fracture networks. For example, fractures represented by one-dimensional fracture elements cannot simulate Mode 2A convection.

Finally, providing a short note on spatial dimensionality is useful. Many studies have been undertaken in two-dimensions but natural systems are three-dimensional where geometry becomes critical. Simply put, in general, all things equivalent, it is easier to get convection in a three-dimensional than a two-dimensional domain. Thus, three-dimensional systems have lower critical *Ra*_{c} than their two-dimensional counterparts. This does not always hold true. For example, for large height/width ratios (vertical slots) the propensity for convection is severely reduced. Further discussion of the effects of spatial dimensionality on free convection is provided by Nield and Bejan (2006), Voss and others (2010), and Knorr and others (2016).