Abstract
Density-driven free convection in porous media is highly affected by large-scale heterogeneity, typical of which are low- or high-permeability inclusions imbedded in homogeneous porous media. In this research, we applied the modified Elder problem to investigate the impact of low- or high-permeability 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 high-permeability inclusion has stronger effects than low-permeability 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-/high-permeability inclusion vertically closer to the source zone, and (4) free convection is more susceptible to the low-permeability inclusion horizontally closer to the source zone. For high-permeability 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 density-driven 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 gravity-induced 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 groundwater-surface 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 local-scale 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 high-permeability inclusions imbedded in homogeneous porous media are particularly common representations of large-scale heterogeneity that can significantly influence solute transport [25]. Post and Simmons [26] used sand tank experiments and numerical models to investigate how discrete low-permeability structures affect free convection at the scale of individual lenses. They found two different mechanisms exist: interlayer convection for high-permeability regions and intralayer convection for low-permeability 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 high-permeability regions.
Moreover, modes of free convection in fractured low-permeability 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 low-permeability 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 high-permeability 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 variable-density flow and solute transport in the system were simulated using SEAWAT-2000 [32]. The permeability and location (horizontal and vertical) of the low- or high-permeability 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-/high-permeability inclusions on free convection. The expected results will advance our understanding of unstable density-drive 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 fixed-concentration boundary at the bottom into a no-flux 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 two-dimensional 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/m3, 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 SEAWAT-2000 [32] to simulate variable-density groundwater flow and solute transport in porous media, which combines MODFLOW-2000 and MT3DMS into a single computer program.
The governing equation for saturated variable-density 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; (L2T-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 gravity-driven 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-1T-1) is the dynamic viscosity; (L2T-1) is the molecular diffusion coefficient; and (L2) is the intrinsic permeability. The value of Ra measures the relation between buoyancy-driven transport and diffusion. The values of the parameters in the equation can be found in Table 1. The critical Ra is for Horton-Rogers-Lapwood 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 Racr, ensuring the occurrence of free convection. Note that when considering a low-/high-permeability 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 no-flux 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 High-Permeability Inclusion
A low-/high-permeability () 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 low-permeability inclusion, 10, 100, and 1000 for a high-permeability 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-/high-permeability 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-/high-permeability 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 time-variable 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 low-permeability inclusion (, Figure 3(a)), the homogeneous case (, Figure 3(b)), and the case of high-permeability inclusion (, Figure 3(c)), it can be observed that after three years five separate fingers were developed in the case of high-permeability region with one main finger generated underneath the high-permeability region in the middle, while only four fingers were generated in the case of low-permeability region and the homogeneous case. This indicates that a higher instability existed in the system of high-permeability region than the other two cases. The high-permeability 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 low-permeability 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 high-permeability inclusion () located at three different horizontal positions are plotted in Figure 4. The salt plumes tend to preferentially develop underneath the high-permeability 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 high-permeability 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 high-permeability 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 higher-permeability 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 high-permeability 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 low-permeability regions () located at three horizontal positions are different (Figure 5) from that for the high-permeability regions. By one year, it can be observed that downward solutes transport in all heterogeneous cases tends to be blocked by the low-permeability 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 low-permeability 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-/high-permeability zone affects free convection. It can be concluded that the high-permeability 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-/high-permeability 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-/high-permeability region were plotted in Figure 6. It can be concluded that the high-permeability region accelerates the mass transport into the domain and the low-permeability 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 high-permeability region is placed vertically further from the source zone (Figure 6(c)). This indicates that the high-permeability region vertically closer to the source zone affects the mass transport in terms of TM more significantly.
For the domain with high-permeability 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 high-permeability 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 high-permeability region () affects the mass transport process, TM values are calculated individually when the high-permeability 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 high-permeability 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 high-permeability region is close to the border of the source zone (, still with the range of the source zone). When the high-permeability 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 high-permeability 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-/high-permeability regions are plotted in Figure 8, in which the Sh value is divided by 60 for better presentation. The results suggested that the high-permeability 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 high-permeability region is placed at the source zone center (position (1, 4), Figure 8(a)). The upgradation is more pronounced when the high-permeability region is placed vertically further from the source zone (Figure 8(c)). Conversely, the low-permeability 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 low-permeability region and preferentially passes through the high-permeability region. As a result, the horizontal COG (COG-) shows right-side biases from the middle of the domain under the presence of the low-permeability region (Figure 9(a)) in the domain and left-side biases under the presence of the high-permeability 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 low-permeability 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 high-permeability 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-/high-permeability 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 low-permeability region in the domain (Figure 9(c)) has less effect on the COG- than that of the high-permeability region (Figure 9(d)). In particular, for the domain with high-permeability 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 density-driven 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-/high-permeability 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 density-driven 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 high-permeability inclusion shows stronger effects on the free convectional process than a low-permeability inclusion at the early stage, with significantly unbalanced solute distributions and accelerated transport processes(2)Free convection is more sensitive to the low-/high-permeability 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 low-permeability inclusion horizontally far from the source zone center. For high-permeability 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-/high-permeability 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.