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 50th 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).
Figure 18 – Elder’s original experimental results in a Hele-Shaw Cell (Elder, 1967b). A fluid is heated from below along part of the base of the Hele-Shaw cell. Elder used suspended aluminum particles to visualize the flow field. The photographs show transient natural convection at two different times: a) six convection cells at early time (T=0.025 dimensionless time); and, b) four convection cells at later time (T=0.05 dimensionless time). One dimensionless time unit (T) in Elder’s experiment is 2.4 × 104 seconds. Thus, (a) with dimensionless time of T=0.025 corresponds to t=10 minutes after the start of the experiment and (b) with T=0.05 corresponds to t=20 minutes after the start of the experiment. In Elder’s original Hele-Shaw cell experiment, plate spacing was 4.0 mm, depth of fluid 5 cm and width of the cavity 20 cm. Reprinted from Journal of Fluid Mechanics, 551, J.W. Elder, Transient convection in a porous medium, 609-623, Copyright (1967), with permission from Cambridge University Press.
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).
Figure 19 – Dimensions, boundary conditions and parameter values of the solute variant of the Elder problem as proposed by Voss and Souza (1987).
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 22l+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 22l+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 (S1), double plume (S2) and triple plume (S3) solutions. Johannsen (2003) found the range of Rayleigh numbers over which these different steady convection solutions could come into being.
Figure 21 – The Elder problem results. Three, steady-state, stable convection states are found using different numerical models and different levels of mesh discretization. These are the single plume (S1), double plume (S2) and triple plume (S3) solutions. Reprinted from Water Resources Research, 45(4), M. Van Reeuwijk and others, Insights from a pseudospectral approach to the Elder problem, Copyright (2009), with permission from John Wiley and Sons.
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 S3, and exists when Ra > 172. A double plume solution, S2, exists when 76 < Ra ≤ 172, as noted by Johannsen (2003). Below Ra = 76, there is only one stable steady state (S1). 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 102 m) in which fingering occurs should always be interpreted with extreme caution.