Research Article  Open Access
Impact of Low or HighPermeability Inclusion on Free Convection in a Porous Medium
Abstract
Densitydriven free convection in porous media is highly affected by largescale heterogeneity, typical of which are low or highpermeability inclusions imbedded in homogeneous porous media. In this research, we applied the modified Elder problem to investigate the impact of low or highpermeability inclusions on the migration of a dense, unstable salt plume. Sensitivity analyses were conducted in terms of the permeability contrast, the effective area (the area of the inclusion beneath the source zone), and the distance of the inclusion from the source zone, all of which were found to play a significant role in controlling the total mass flux released from the source into the media. Results show that (1) a highpermeability inclusion has stronger effects than lowpermeability inclusion, due to significantly unbalanced solute distributions caused by accelerated solute transport, (2) the inclusion with a larger effective area has more potential to influence free convection, (3) free convection is more sensitive to the low/highpermeability inclusion vertically closer to the source zone, and (4) free convection is more susceptible to the lowpermeability inclusion horizontally closer to the source zone. For highpermeability inclusions, the inclusion horizontally closer to the source zone influences the transport process more significantly at the early stage, and conversely, the inclusion far from the source zone has a later impact. The results obtained could offer significant implications for understanding unstable densitydriven flow and solute transport in porous media with structured heterogeneity.
1. Introduction
In groundwater systems, the fluid density gradient may vary as a function of concentration, temperature, and/or pressure of the fluid, which is of particular importance in the migration of solutes. The density gradient caused by a denser fluid overlying a light one may lead to the onset of the gravityinduced instabilities. Firstly, a boundary layer develops at the source zone through diffusion and grows to a certain critical thickness and then breaks into fingers moving downward rapidly, generating the instability (free convection) [1–3]. Over the years, numerical, laboratory, and field experiments have analyzed numerous examples of naturally occurring free convection, including seawater intrusion in coastal aquifers, leachate infiltration from waste disposal, salt lakes and saline disposal basins, and interaction between groundwatersurface water [e.g., [4–12]]. Free convection significantly enhances the hydrodynamic mixing of denser solute with its ambient groundwater over lager spatial scales within shorter time scales compared with diffusion alone.
Heterogeneity in porous media can significantly affect the transport processes of free convection. In addition to the common assumption of homogeneous and isotropic hydrogeological settings in the studies of free convection processes (for exhaustive reviews on this topic, see [13–16]), heterogeneity has also been frequently investigated in previous studies. [8] conducted an experimental investigation of free convection in homogeneous, layered, and lenticular porous media to study the roles played by the heterogeneity. They noted that localscale nonhomogeneity of the flow results in a complex and highly dispersed plume. [17] found that in flow through natural porous media, heterogeneity can cause perturbation or interfacial disturbances over many spatial scales (from slight differences in pore geometry to larger heterogeneities on the scale of the dense plume source). [18] improved our understanding of instability phenomena that increased correlation length scales and increased variance of the permeability field could promote stability. The experimental and theoretical studies mentioned above have proved that the passive solute transport in natural porous media is controlled by the spatial variations in heterogeneous fields. Another notable result in the field of free convection was the quantitative indicators developed by [19], which include solute fluxes, solute presence, center of gravity, and finger penetration depth, providing visual inspection to assess numerical outputs (e.g., [20–23]). [15] carried out a numerical investigation of free convection and considered various patterns of heterogeneity ranging from periodically structured to random permeability fields. The results suggested that heterogeneity triggers the onset of instabilities and also tend to dissipate them at later times by enhancing dispersive mixing. They also found that structured heterogeneity allows instabilities to propagate at a modest combination of fracture aperture and separation distances, while disordered heterogeneity tends to dissipate convection through dispersive mixing. [24] evaluated free thermal convection in randomly heterogeneous media.
In heterogeneous porous media, low or highpermeability inclusions imbedded in homogeneous porous media are particularly common representations of largescale heterogeneity that can significantly influence solute transport [25]. Post and Simmons [26] used sand tank experiments and numerical models to investigate how discrete lowpermeability structures affect free convection at the scale of individual lenses. They found two different mechanisms exist: interlayer convection for highpermeability regions and intralayer convection for lowpermeability regions. The result demonstrates that the solute transport driven by a source zone at the top of the domain is a complicated function of the geometry of the permeability structure and contrast between low and highpermeability regions.
Moreover, modes of free convection in fractured lowpermeability media have been explored by [27]. Their results showed that free convection may be common in thick shale sequences, and analyses neglecting both inter and intrafracture convection modes may significantly underestimate the likelihood of the occurrence of free convection. [28] studied fracture networks containing continuous, discontinuous, orthogonal, and/or inclined discrete fractures embedded in a lowpermeability rock matrix and found that fracture networks can inhibit or promote convection depending on the fracture network geometry. [29] investigated combined inter and intrafracture free convection in fractured porous rock and suggested that it is challenging to predict the convective pattern under such conditions.
In this paper, we focus on investigating the effects of the distribution of low or highpermeability inclusions on free convection. We introduce structured heterogeneity in porous media in a modified setting of the Elder problem [30, 31], which is a commonly used example of free convection. The variabledensity flow and solute transport in the system were simulated using SEAWAT2000 [32]. The permeability and location (horizontal and vertical) of the low or highpermeability inclusions were set to be variable. Their effects on free convection were quantitatively assessed using three assessment indicators, including the total solute mass (TM), Sherwood number (Sh), and solute center of gravity (COG). The main objective of the current study was to test the sensitivity of the location of the low/highpermeability inclusions on free convection. The expected results will advance our understanding of unstable densitydrive flow and solute transport processes under the effect of structured heterogeneity.
2. Methodology
2.1. Modified Elder Problem
The Elder problem is a classic example of free convection phenomena [33, 34]. It was first studied by [30, 31] for thermal convection (referred as “original” Elder problem) and developed later by [35] into a solute transport problem (referred as “classic” Elder problem). The Elder problem is chosen here for two reasons: (1) it is a classic example for free convection in which the flow is purely driven by density gradients and (2) the Elder problem, especially under the homogeneous condition, has been extensively investigated and discussed in both experimental and numerical simulation researches.
The Elder problem used in this study (Figure 1) is modified from the one used in [21], by transforming the fixedconcentration boundary at the bottom into a noflux boundary while keeping the domain geometry and boundary conditions unchanged. This transformation makes the solute transport process more realistic under natural conditions. Initially, the domain is filled with freshwater with a hydrostatic head distribution. The twodimensional domain has a length of 600 m and a depth of 150 m. The domain boundaries are impermeable except that constant heads are prescribed at the top left and right corners as two outlets for fluid and solute. A fixed concentration boundary condition ( [20, 23]) is imposed at the source zone. The density of the source zone () is 1200 kg/m^{3}, which corresponds to a salt concentration of . As abundant solutes diffused from the source zone accumulate at the top boundary, flow is generated by density gradients in the model domain.
2.2. Numerical Simulations
2.2.1. Simulation Tool
As indicated above, we employ SEAWAT2000 [32] to simulate variabledensity groundwater flow and solute transport in porous media, which combines MODFLOW2000 and MT3DMS into a single computer program.
The governing equation for saturated variabledensity groundwater flow in terms of freshwater head is expressed by [36] where (L) is the vertical coordinate directed upward; (L) is the equivalent freshwater head; (LT^{1}) is the equivalent freshwater hydraulic conductivity; (ML^{3}), (ML^{3}), and (ML^{3}) are the densities of fluid, freshwater, and source/sink, respectively; () is effective porosity; (L^{1}) is the equivalent freshwater storage coefficient; (T) is the time; and (T^{1}) is the flow rate per unit volume of aquifer of the source/sink.
The solute transport model involving advection, molecular diffusion, and mechanical dispersion is given as where (ML^{3}) is dissolved concentration; (L^{2}T^{1}) is the hydrodynamic dispersion coefficient tensor; and (LT^{1}) is the pore water velocity vector.
An empirical linear relationship is present between salt concentration and density (e.g., [37, 38]): where () has a dimensionless value of 0.7143 [32].
2.2.2. Instability Descriptor
Stratification formed by a denser fluid overlying a lighter fluid develops high a density gradient and causes gravitydriven instability. The instability of the Elder problem can be quantified by the dimensionless Rayleigh number [15]: where is the dimensionless density contrast coefficient; (L) is the height of the model domain; (LT^{2}) is the gravitational acceleration; (ML^{1}T^{1}) is the dynamic viscosity; (L^{2}T^{1}) is the molecular diffusion coefficient; and (L^{2}) is the intrinsic permeability. The value of Ra measures the relation between buoyancydriven transport and diffusion. The values of the parameters in the equation can be found in Table 1. The critical Ra is for HortonRogersLapwood problem with infinitely extended horizontal porous layer and constant temperature at upper and lower boundaries, but is zero for the classic Elder problem [23, 39] since the fixed concentration boundary at both upper and lower boundaries decides the inevitable concentration gradient and thus free convection. The calculated Rayleigh number of our modified Elder problem approaches 400, which is far greater than the critical value Ra_{cr}, ensuring the occurrence of free convection. Note that when considering a low/highpermeability inclusion in the domain, the value of Ra alters accordingly.

2.2.3. Spatial and Temporal Discretization
Previous studies have proved that free convection processes in the Elder problem are largely influenced by grid discretization [19, 40]. [18] showed that too coarse grid discretization may result in large errors and numerical dispersion, which can reshape the salt plumes and transform the developments of unstable fingers after their generation. As suggested by [41], (horizontal) and (vertical) are adopted in this study to avoid the influences of grid discretization, resulting in a total of rectangular cells.
As the bottom boundary condition of the modified Elder problem is transformed to be a noflux boundary, the system cannot reach its steady state until the entire domain is filled with saltwater with a concentration equal to that of the source zone. Therefore, the simulation time is extended to 18,250 days (50 yrs), such that a significant amount of solute mass could enter into the domain.
2.2.4. Low or HighPermeability Inclusion
A low/highpermeability () inclusion is imbedded in the modeled porous media of constant permeability (). The permeability contrast is set to be 0.01 and 0.1 for a lowpermeability inclusion, 10, 100, and 1000 for a highpermeability inclusion, and 1 for the homogeneous porous media as a reference case.
Due to the symmetry of the model domain, only the left half of the domain was subjected to the simulations (Figure 2). Sixteen zones were selected for the low/highpermeability inclusion, each with a rectangular area of . Their locations, labeled by their horizontal and vertical positions (Figure 2), were selected to assess the behavior of the free convection in response to the changing location of the low (high)permeability inclusion. Note that the leftmost two inclusions are partially overlapped. Preliminary simulations showed that the low (high)permeability inclusions located close to the top boundary affect the free convection mostly. Therefore, emphasis was more put on the low (high)permeability inclusions situated in regions from (1, 1) to (1, 4) (Figure 2). For cases of , 100, and 1000, four positions of inclusions are considered (i.e., from (1, 1) to (1, 4)), while for cases of and 10, 16 positons are adopted (i.e., from (1, 1) to (4, 4)). Table 2 lists the cases with low (high)permeability inclusions. By considering a homogeneous case, a total of 45 cases are simulated.

2.3. Assessment Indicators
To quantify the behavior of free convection in response to the permeability/location of the low/highpermeability inclusions, three assessment indicators [19, 21–23] are used in this study.
2.3.1. Total Solute Mass (TM)
TM is a parameter to quantify timevariable solute mass (normalized) present in the model domain, calculated as [21] where is the normalized concentration ranging within [0, 1] at (). The value of TM varies from 0 when the domain has no solute mass to 1 when the domain is fully filled with saltwater.
2.3.2. Sherwood Number (Sh)
Sh is defined as the ratio of actual mass transfer due to free convection during the transient state to the rate of mass transfer due to diffusion, calculated as [14] where is the solute mass flux across the top boundary, and are the width and length of the source zone, and () is the maximum concentration difference between freshwater and saltwater. A stable system with only diffusive transport has . Conversely, an unstable system with convection and diffusion, such as the classic Elder problem, has .
2.3.3. Center of Gravity (COG) of the Plume
COG describes the horizontal and vertical transient gravity center of solute plume present in the domain. It is one of the most important stochastic characteristics and can be calculated through spatial moments: where the summation is executed for all cells of the model grid; is the normalized concentration of cell (); and are coordinates of the cell; and and are integer power exponents that define the moment order. Subsequently, the gravity center () is expressed by
COG can indicate the movement of the salt plumes in the heterogeneous porous media in both vertical and horizontal direction.
3. Results and Discussion
3.1. Plume Development
Comparing the concentration distributions for the case of lowpermeability inclusion (, Figure 3(a)), the homogeneous case (, Figure 3(b)), and the case of highpermeability inclusion (, Figure 3(c)), it can be observed that after three years five separate fingers were developed in the case of highpermeability region with one main finger generated underneath the highpermeability region in the middle, while only four fingers were generated in the case of lowpermeability region and the homogeneous case. This indicates that a higher instability existed in the system of highpermeability region than the other two cases. The highpermeability region aggravates the imbalance by accelerating its transport in the local region, resulting in the deeper and wider penetration of the salt plumes. Conversely, the lowpermeability region slows down the process of descending and spreading, such that the solute distributes more uniformly. Conclusion can be made that, in terms of instability, a higher permeability region in the domain has a pronounced impact on free convection, which would result in more solute mass entering into the domain at the same loading time. Hence, emphasis for discussion will be put on the cases with higher permeability region () afterwards.
The concentration distributions in the domain with the highpermeability inclusion () located at three different horizontal positions are plotted in Figure 4. The salt plumes tend to preferentially develop underneath the highpermeability region, consistent with our aforementioned conclusion that a higher permeability region can result in a faster flow speed and a less uniform concentration distribution. After one year, an additional finger is generated when the highpermeability region is located within the source zone (Figure 4(c)). By three years, except case (1, 1) with a less effective area, the main finger developed from the highpermeability region starts to reach the domain bottom, earlier than that in the homogeneous case. Surprisingly, at approximately 15 years, all fingers have merged into a single dense plume. This may be because the area of the higherpermeability inclusion is relatively small in comparison to the total domain area. The development of the salt plumes under all cases experiences such a process, i.e., individual fingers developed firstly and merged into one large plume later. However, diversity existed in the solute migration paths that depended on the position of the highpermeability region. Eventually, as the density gradient dissipates, diffusion becomes the primary transport mechanism, indicating that the system state was gradually transformed from instability to stability.
The concentration distributions for the domain with lowpermeability regions () located at three horizontal positions are different (Figure 5) from that for the highpermeability regions. By one year, it can be observed that downward solutes transport in all heterogeneous cases tends to be blocked by the lowpermeability region, such that only one finger was developed in two of the cases (Figures 5(a) and 5(b)). Three years after, the distributions show that part of salt plumes immigrate bypass the lowpermeability regions rather than passing through them. Fifty years after, the systems approached their steady state with the plumes symmetrically distributed.
Insights into the plume development have provided qualitative explanations for the ways how the presence of a spatially varied low/highpermeability zone affects free convection. It can be concluded that the highpermeability region () in the domain has more potential to influence the transport process of free convection in time and space (e.g., fingering speed, number, and pattern). Quantitative evaluation of the behavior of the system in response to the presence of a low/highpermeability region requires further investigation using the assessment indicators.
3.2. Quantitative Assessment
3.2.1. Effects on Total Solute Mass (TM)
The breakthrough curves for TM of the model domain with different configurations of the low/highpermeability region were plotted in Figure 6. It can be concluded that the highpermeability region accelerates the mass transport into the domain and the lowpermeability region slows down that transport processes, evidenced by TM curves overtopping (Figure 6(a)) and underneath (Figure 6(b)) the one for homogeneous case, respectively. The TM curves almost coincide with each other when the highpermeability region is placed vertically further from the source zone (Figure 6(c)). This indicates that the highpermeability region vertically closer to the source zone affects the mass transport in terms of TM more significantly.
For the domain with highpermeability regions close to the source zone (Figure 6(a)), linear increasing trends of the TM can be identified at an early stage (<5 yr) while the big main fingers develop under high concentration gradients (Figure 4). Additionally, the transient TM is higher when the highpermeability region is closer to the source zone (position (1, 4), Figure 6(a)) at the early stage. The TM increasing speeds slow down gradually when salt is filling the domain after the stage of finger developing.
To test how the horizontal distance between the source zone center and the highpermeability region () affects the mass transport process, TM values are calculated individually when the highpermeability region is horizontally moved from the source zone center (, Figure 7) to the side of the source zone (, Figure 7). At the early stage, e.g., , the TM shows a monotonic decreasing function, i.e., mass transports faster when the highpermeability region is close to the source zone due to easier finger generation. However, this trend does not exist anymore at the later stage, within which the maximum TM is found when the highpermeability region is close to the border of the source zone (, still with the range of the source zone). When the highpermeability region is placed further sideward out of the range of the source zone (), the acceleration effect on mass transport process becomes less pronounced. In summary, the distance between the highpermeability region and the source zone significantly determines their effect on the transport process in terms of total mass. The trends differ between different stages of the transport process. Initially, the location of the inclusion dominates the solute loading process, while the total mass in the domain controls the speed of solute entering in to the domain at larger times. The time dependence of the TM on the location of inclusion is similar to the effect of dual domain mass transfer found in our previous study [20].
3.2.2. Effects on Sherwood Number (Sh)
The calculated Sh values for the domain with different configurations of the low/highpermeability regions are plotted in Figure 8, in which the Sh value is divided by 60 for better presentation. The results suggested that the highpermeability regions have upgraded the convectional mass transport, proved by higher Sh values than the homogeneous case (Figure 8(a)). At the early stage, this upgrading effect is more pronounced when the highpermeability region is placed at the source zone center (position (1, 4), Figure 8(a)). The upgradation is more pronounced when the highpermeability region is placed vertically further from the source zone (Figure 8(c)). Conversely, the lowpermeability regions tend to degrade the convectional mass transport with lower Sh numbers than the homogenous case (Figure 8(b)). However, these upgradation or degradation of convectional mass transport becomes less significant at the later stage with the Sh values of all cases converging to a single value.
3.2.3. Effects on the Center of Gravity (COG)
Nonuniform concentration distributions and fingering processes can certainly result in horizontal and vertical bias of the COG. Generally, solutes transport preferentially takes the path through more permeable area, i.e., salt plumes bypass the lowpermeability region and preferentially passes through the highpermeability region. As a result, the horizontal COG (COG) shows rightside biases from the middle of the domain under the presence of the lowpermeability region (Figure 9(a)) in the domain and leftside biases under the presence of the highpermeability region (Figure 9(b)). Note that the COG shows no difference from the homogenous case when the low (high)permeability regions are placed in the middle of the domain (position (1, 4), Figure 9) due to the symmetric distribution of the salt plumes. For the domain with lowpermeability regions (Figure 9(a)), fingers preferentially develop on the right side of the domain, with the bias of the COG reaching their maximums when the main fingers touch the domain bottom. Later, the plumes start to develop horizontally with the bias of COG being continuously reduced. However, double peaks of the bias of COG can be observed for the cases of highpermeability regions (Figure 9(b)), with the first one reached when the main finger touches the domain bottom and the second one reached due to a large amount of solute accumulation and diffusion on the left side of the domain. Additionally, the low/highpermeability region placed the farthest from the source zone (but still within the source zone, position (1, 2), Figure 9) causes the most significant bias of COG. At the late stage, the COG moves continuously towards the domain center as secondary fingers go reach the bottom and salt plumes spread by diffusion to the entire space of the domain.
The vertical COG (COG) for all cases shows that the stage of fast descend due to the penetration of the fingers, followed by a stable stage during which the COG slightly rebounds since the fingers touch the domain bottom. The location of the lowpermeability region in the domain (Figure 9(c)) has less effect on the COG than that of the highpermeability region (Figure 9(d)). In particular, for the domain with highpermeability region in the domain, deeper COG can be reached when the region is closer to the source (position (1, 4), Figure 9(d)) due to less lateral transport of the solute at the early stage.
4. Summary and Conclusion
In this study, we focused on investigating the influence of permeability inclusions on densitydriven free convection. We introduced the structured heterogeneity in porous media for a modified setting of the Elder problem. The effects of the permeability and the location (horizontal and vertical) of the low/highpermeability inclusion on the free convection process were assessed using three assessment indicators: total solute mass (TM), Sherwood number (Sh), and solute center of gravity (COG). Despite the simplified model assumed, the obtained numerical results have important implications for understanding the unstable densitydriven flow and solute transport processes in the presence of structured heterogeneity. The key findings of the study are as follows: (1)In comparison to homogeneous systems, heterogeneity in permeability makes the system more unstable. A highpermeability inclusion shows stronger effects on the free convectional process than a lowpermeability inclusion at the early stage, with significantly unbalanced solute distributions and accelerated transport processes(2)Free convection is more sensitive to the low/highpermeability inclusion vertically closer to the source zone, revealed in the mostly influenced breakthrough curves for total solute mass and Sherwood number(3)Free convection is more susceptible to the lowpermeability inclusion horizontally far from the source zone center. For highpermeability inclusions, by contrast, the inclusion horizontally closer to the source zone center influences the convectional transport process more significantly at the early stage, while the inclusion far from the source zone center has more significant influences at the late stage
The results obtained in the current study are based on the Elder problem, which neglects the boundary conditions in real situations. For example, in a coastal setting, one of the vertical boundaries should be the coastal boundary. Under such conditions, a constant head inland boundary with a constant salt concentration is usually defined, resulting in seawater intrusion in a coastal aquifer. The consideration of real situations incorporated with the effect of a low/highpermeability inclusion should be included in future studies.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
C. Lu acknowledges the financial support from the National Key Research Project (2018YFC0407200), the National Natural Science Foundation of China (51679067 and 51879088), the Fundamental Research Funds for the Central Universities (2018B42814), and the “111” Project (B17015). Y. Xie acknowledges the financial support from the National Natural Science Foundation of China (41702243). All the results/data have been included in the paper.
References
 R. A. Wooding, “Free convection of fluid in a vertical tube filled with porous material,” Journal of Fluid Mechanics, vol. 13, no. 1, pp. 129–144, 1962. View at: Publisher Site  Google Scholar
 R. A. Wooding, S. W. Tyler, and I. White, “Convection in groundwater below an evaporating Salt Lake: 1. Onset of instability,” Water Resources Research, vol. 33, no. 6, pp. 1199–1217, 1997. View at: Publisher Site  Google Scholar
 R. A. Wooding, S. W. Tyler, I. White, and P. A. Anderson, “Convection in groundwater below an evaporating Salt Lake: 2. Evolution of fingers or plumes,” Water Resources Research, vol. 33, no. 6, pp. 1219–1228, 1997. View at: Publisher Site  Google Scholar
 E. O. Frind, “Simulation of longterm transient densitydependent transport in groundwater,” Advances in Water Resources, vol. 5, no. 2, pp. 73–88, 1982. View at: Publisher Site  Google Scholar
 P. S. Huyakorn, P. F. Andersen, J. W. Mercer, and H. O. White Jr., “Saltwater intrusion in aquifers: development and testing of a threedimensional finite element model,” Water Resources Research, vol. 23, no. 2, pp. 293–312, 1987. View at: Publisher Site  Google Scholar
 H. Kooi, J. Groen, and A. Leijnse, “Modes of seawater intrusion during transgressions,” Water Resources Research, vol. 36, no. 12, pp. 3581–3589, 2000. View at: Publisher Site  Google Scholar
 H. A. Michael, K. C. Scott, M. Koneshloo, X. Yu, M. R. Khan, and K. Li, “Geologic influence on groundwater salinity drives large seawater circulation through the continental shelf,” Geophysical Research Letters, vol. 43, no. 20, pp. 10,782–10,791, 2016. View at: Publisher Site  Google Scholar
 R. A. Schincariol and F. W. Schwartz, “An experimental investigation of variable density flow and mixing in homogeneous and heterogeneous media,” Water Resources Research, vol. 26, no. 10, pp. 2317–2329, 1990. View at: Publisher Site  Google Scholar
 R. A. Schincariol, F. W. Schwartz, and C. A. Mendoza, “On the generation of instabilities in variable density flow,” Water Resources Research, vol. 30, no. 4, pp. 913–927, 1994. View at: Publisher Site  Google Scholar
 C. T. Simmons and K. A. Narayan, “Mixed convection processes below a saline disposal basin,” Journal of Hydrology, vol. 194, no. 1–4, pp. 263–285, 1997. View at: Publisher Site  Google Scholar
 J. D. Stevens, J. M. Sharp Jr., C. T. Simmons, and T. R. Fenstemaker, “Evidence of free convection in groundwater: fieldbased measurements beneath windtidal flats,” Journal of Hydrology, vol. 375, no. 34, pp. 394–409, 2009. View at: Publisher Site  Google Scholar
 R. L. Van Dam, C. T. Simmons, D. W. Hyndman, and W. W. Wood, “Natural free convection in porous media: first field documentation in groundwater,” Geophysical Research Letters, vol. 36, no. 11, 2009. View at: Publisher Site  Google Scholar
 B. Gebhart, Y. Jaluria, R. L. Mahajan, and B. Sammakia, BuoyancyInduced Flows and Transport, Hemisphere Publishing Corp., Washingtong DC., 1988.
 D. A. Nield and A. Bejan, Convection in Porous Media, SpringerVerlag, New York, 2013. View at: Publisher Site
 C. T. Simmons, T. R. Fenstemaker, and J. M. Sharp Jr., “Variabledensity groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges,” Journal of Contaminant Hydrology, vol. 52, no. 14, pp. 245–275, 2001. View at: Publisher Site  Google Scholar
 C. T. Simmons, “Variable density groundwater flow: from current challenges to future possibilities,” Hydrogeology Journal, vol. 13, no. 1, pp. 116–119, 2005. View at: Publisher Site  Google Scholar
 D. E. Moissis and M. F. Wheeler, “Effect of the structure of the porous medium on unstable miscible displacement,” in Dynamics of Fluids in Hierachical Porous Medida, J. H. Cushman, Ed., pp. 243–271, Academic, San Diego, California, 1990. View at: Google Scholar
 R. A. Schincariol, F. W. Schwartz, and C. A. Mendoza, “Instabilities in variable density flows: stability and sensitivity analyses for homogeneous and heterogeneous media,” Water Resources Research, vol. 33, no. 1, pp. 31–41, 1997. View at: Publisher Site  Google Scholar
 A. Prasad and C. T. Simmons, “Unstable densitydriven flow in heterogeneous porous media: A stochastic study of theElder[1967b] “short heater” problem,” Water Resources Research, vol. 39, no. 1, pp. SBH 41–SBH 421, 2003. View at: Publisher Site  Google Scholar
 C. Lu, L. Shi, Y. Chen, Y. Xie, and C. T. Simmons, “Impact of kinetic mass transfer on free convection in a porous medium,” Water Resources Research, vol. 52, no. 5, pp. 3637–3653, 2016. View at: Publisher Site  Google Scholar
 V. E. A. Post and H. Prommer, “Multicomponent reactive transport simulation of the Elder problem: effects of chemical reactions on salt plume development,” Water Resources Research, vol. 43, no. 10, 2007. View at: Publisher Site  Google Scholar
 A. Prasad and C. T. Simmons, “Using quantitative indicators to evaluate results from variabledensity groundwater flow models,” Hydrogeology Journal, vol. 13, no. 56, pp. 905–914, 2005. View at: Publisher Site  Google Scholar
 Y. Xie, C. T. Simmons, A. D. Werner, and J. D. Ward, “Effect of transient solute loading on free convection in porous media,” Water Resources Research, vol. 46, no. 11, 2010. View at: Publisher Site  Google Scholar
 V. T. Nguyen, T. Graf, and C. R. Guevara Morel, “Free thermal convection in heterogeneous porous media,” Geothermics, vol. 64, pp. 152–162, 2016. View at: Publisher Site  Google Scholar
 J. Luo and O. A. Cirpka, “How well do mean breakthrough curves predict mixingcontrolled reactive transport?” Water Resources Research, vol. 47, no. 2, 2011. View at: Publisher Site  Google Scholar
 V. E. A. Post and C. T. Simmons, “Free convective controls on sequestration of salts into lowpermeability strata: insights from sand tank laboratory experiments and numerical modelling,” Hydrogeology Journal, vol. 18, no. 1, pp. 39–54, 2010. View at: Publisher Site  Google Scholar
 C. T. Simmons, J. M. Sharp Jr., and D. A. Nield, “Modes of free convection in fractured lowpermeability media,” Water Resources Research, vol. 44, no. 3, 2008. View at: Publisher Site  Google Scholar
 K. Vujević, T. Graf, C. T. Simmons, and A. D. Werner, “Impact of fracture network geometry on free convective flow patterns,” Advances in Water Resources, vol. 71, pp. 65–80, 2014. View at: Publisher Site  Google Scholar
 K. Vujevic and T. Graf, “Combined inter and intrafracture free convection in fracture networks embedded in a lowpermeability matrix,” Advances in Water Resources, vol. 84, pp. 52–63, 2015. View at: Publisher Site  Google Scholar
 J. W. Elder, “Transient convection in a porous medium,” Journal of Fluid Mechanics, vol. 27, no. 3, pp. 609–623, 1967. View at: Publisher Site  Google Scholar
 J. W. Elder, “Steady free convection in a porous medium heated from below,” Journal of Fluid Mechanics, vol. 27, no. 1, pp. 29–48, 1967. View at: Publisher Site  Google Scholar
 C. D. Langevin, W. B. Shoemaker, and W. Guo, “Modflow2000, the US geological survey modular groundwater modeldocumentation of the SEAWAT2000 version with the variabledensity flow process (VDF) and the integrated MT3DMS transport process (IMT),” U.S. Geological Survey OpenFile Report, vol. 03426, p. 43, 2003. View at: Google Scholar
 J. W. Elder, C. T. Simmons, H. J. Diersch, P. Frolkovič, E. Holzbecher, and K. Johannsen, “The Elder problem,” Fluids, vol. 2, no. 1, p. 11, 2017. View at: Publisher Site  Google Scholar
 C. T. Simmons and J. W. Elder, “The Elder problem,” Groundwater, vol. 55, no. 6, pp. 926–930, 2017. View at: Publisher Site  Google Scholar
 C. I. Voss and W. R. Souza, “Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwatersaltwater transition zone,” Water Resources Research, vol. 23, no. 10, pp. 1851–1866, 1987. View at: Publisher Site  Google Scholar
 C. D. Langevin and W. Guo, “MODFLOW/MT3DMSbased simulation of variabledensity ground water flow and transport,” Groundwater, vol. 44, no. 3, pp. 339–351, 2006. View at: Publisher Site  Google Scholar
 M. J. Simpson and T. P. Clement, “Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of densitydependent groundwater flow models,” Advances in Water Resources, vol. 26, no. 1, pp. 17–31, 2003. View at: Publisher Site  Google Scholar
 H. Zhang and F. W. Schwartz, “Multispecies contaminant plumes in variable density flow systems,” Water Resources Research, vol. 31, no. 4, pp. 837–847, 1995. View at: Publisher Site  Google Scholar
 M. van Reeuwijk, S. A. Mathias, C. T. Simmons, and J. D. Ward, “Insights from a pseudospectral approach to the Elder problem,” Water Resources Research, vol. 45, no. 4, 2009. View at: Publisher Site  Google Scholar
 H. J. G. Diersch and O. Kolditz, “Variabledensity flow and transport in porous media: approaches and challenges,” Advances in Water Resources, vol. 25, no. 8–12, pp. 899–944, 2002. View at: Publisher Site  Google Scholar
 A. AlMaktoumi, D. A. Lockington, and R. E. Volker, “SEAWAT 2000: modelling unstable flow and sensitivity to discretization levels and numerical schemes,” Hydrogeology Journal, vol. 15, no. 6, pp. 1119–1129, 2007. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Min Yan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.