# 6.2 The Elder Problem

Because the governing equations are non-linear and coupled, analytical solutions are intractable and only exist for the simplest arrangements of convection cells. This means that numerical codes are not easily verified, i.e., the testing of their accuracy by comparing the numerical result to an analytical solution is generally not possible. In that case, alternative strategies are required, such as internal consistency testing (for example mass balance checks) and external consistency testing (Konikow et al., 1996). The latter involves comparing the results of different numerical codes against one another, which is known as benchmarking. A number of model benchmarks have been developed and applied over the years to test variable-density flow models. Information on these model benchmarks and test cases, including successes and limitations, can be found in the references noted in the following discussion and review papers (Diersch and Kolditz, 2002; Voss et al., 2010 and Simmons et al., 2010).

One of the most intensely studied and well-known benchmark problems for variable-density simulators is the Elder problem. It is based on pioneering work by John W. Elder in the Cavendish Laboratory at Cambridge. It was published in the Journal of Fluid Mechanics in 1967 (Elder, 1967b, 1967a), and two papers marking the 50^{th} anniversary of Elder’s work were published in 2017 (Elder et al., 2017; Simmons and Elder, 2017). What is remarkable is that he did not realize his work at Cambridge had gone on to become the basis for the benchmark problem that became famous in hydrogeology. Elder’s original study was based on a thermal experiment conducted in a Hele-Shaw Cell. A Hele-Shaw cell consists of two closely spaced parallel plates, between which a fluid, or multiple fluids, are injected. Because the flow in the cell obeys Darcy’s law, Hele-Shaw cells can be used to study flow in porous media (having an intrinsic permeability equal to where *a* is the distance between the plates). Dispersive transport, however, cannot be studied (NEA, 1990).

In Elder’s experiment, the fluid was heated from below across part of its base (a “short heater”). The fluid flow was visualized using suspended aluminum particles and the photographs of the results are shown in Figure 18 (Elder, 1967b). Six convection cells occur at early time (Figure 18a) and form four larger cells at later time (Figure 18b).

The solute analog was developed by authors including Diersch (1981), Voss and Souza (1987), and Holzbecher (1991). Figure 19 shows the boundary conditions for the solute analog of the Elder problem. Rather than heating from below, the solute analog has solute entering from above. The Rayleigh number is *Ra = *400. It was considered part of the OECD (Organisation for Economic Co-operation and Development) Hydrocoin Project which was designed to test numerical simulators as part of the modeling required for nuclear waste repositories (NEA, 1992).

The Hydrocoin Project made it clear that as different numerical codes were used to simulate the Elder problem, their results were generally similar, but there were differences in the details. That is, even though the codes roughly predicted the same fingering pattern, the solute concentration fields were not exactly the same. Moreover, the project concluded that the results were sensitive to the grid discretization (NEA, 1992, 1990). Later on, both points became a topic of intensive research, which revealed that the differences between the models occurred to the differences in the nature and number of cells and whether there was a central downwelling or upwelling (Diersch and Kolditz, 2002).

Figure 20 gives an example of varying results as a function of the level of numerical grid resolution. The results were calculated using SEAWAT (Langevin and Guo, 2006) based on the grid refinement approach by Frolkovič and De Schepper (2000). The grid refinement, *l*, determines the number of nodes according to 2^{2l+1}. As the grid is refined, the nature and number of convection cells changes, and there can be upwelling or downwelling below the center of the “short heater”. This variation in the flow direction was summarized by Diersch and Kolditz (2002). They found that different authors reported different results depending on the numerical code and the level of mesh discretization.

**Figure ****20** – The Elder problem results (*Ra*=400) as grid level, *l*, is varied. Small *l* is a coarser grid and larger *l* is finer grid. The number of cells is 2^{2l+1}. Lower rows represent later time.

Frolkovič and De Schepper (2000) noted that their models evolved to one of three plume configurations as they were run long enough to reach steady solutions. These are shown in Figure 21 and are called the single plume (S_{1}), double plume (S_{2}) and triple plume (S_{3}) solutions. Johannsen (2003) found the range of Rayleigh numbers over which these different steady convection solutions could come into being.

The question that arose was if these different steady-state convection solutions were real or numerical artifacts? To resolve the issue, Van Reeuwijk and others (2009) used a pseudospectral method to solve the Elder problem. This method eliminates all sources of numerical discrepancy, to obtain uncontaminated insights. The use of a pseudospectral method avoids all truncation errors associated with differentiation. The Van Reeuwijk and others (2009) results confirmed the earlier results from Johannsen (2003) and Frolkovič and De Schepper (2000): at *Ra* = 400 three stable steady-state solutions co-exist. Hence, the ambiguities are physical rather than numerical. This regime is called the triple plume solution, denoted S_{3}, and exists when *Ra* > 172. A double plume solution, S_{2}, exists when 76 < *Ra* ≤ 172, as noted by Johannsen (2003). Below *Ra* = 76, there is only one stable steady state (S_{1}). Van Reeuwijk and others (2009) therefore suggested the use of a low Rayleigh number solution for benchmark studies. This particular variant of the Elder problem, which has *Ra* = 60, is called the Low Rayleigh Number Elder problem.

Grid convergence is used to assess the performance of groundwater models of constant-density systems where the numerical errors tend to get smaller as the grid is refined. The studies by Van Reeuwijk and others (2009), Johannsen (2003) and Frolkovič and De Schepper (2000) led to the insight that the usual notion of grid convergence does not apply to the Elder problem, nor to numerical models of free convection models more generally, because it is impossible to expect to find a single answer when the physical system being simulated can reach more than one state under identical conditions, that is, there are multiple valid solutions.

The studies of the Elder and other benchmark problems have shown that there can be significant variability in results for the same problem when using different numerical software codes. The results of numerical simulations of free convection are sensitive to solution schemes used in the model. One critical aspect is the generation of heterogeneities that control the onset of convective fingering. Small perturbations of the concentration or temperature are required to trigger the unstable behavior, and in a numerical code these are caused by round-off error unless they are somehow prescribed by the modeler. Knowing how to prescribe perturbations is not trivial and the result of the simulation is inevitably dependent to at least some extent on the method selected for introducing perturbations.

In addition to the aforementioned issues associated with discretization, the physical dimensions of the boundary layer and the plume dimensions limit the maximum model cell dimensions (Kooi et al., 2000b). Fingers form once a boundary layer reaches a critical thickness, hence the vertical discretization of a numerical model must be such that the cell height, *∆**z*, is less than *δ*_{cr}, where *δ*_{cr} is the critical boundary layer thickness (Section 5.3). With *∆**z* > *δ*_{cr} a boundary layer is created that is numerically unstable, hence the onset of fingering is simulated unrealistically early. Moreover, the cell width *∆**x* must be smaller than the instabilities that form, the size of which depends on the critical boundary layer thickness. For example, for a boundary layer that is compressed by an upward flux (Section 5.3), (Kooi et al., 2000b) recommended that *∆**x* be smaller than the critical instability wavelength, hence they presented Equation 37 in contrast to Equation 29.

(37) |

The fact that *k* is in the denominator can mean that* ∆**x* can become impractically small for permeable layers (Post and Kooi, 2003). For this reason, the results of regional-scale models (with model grid cell dimensions typically on the order of 10 to 10^{2} m) in which fingering occurs should always be interpreted with extreme caution.