Convective Flow in an Aquifer Layer

Here, we investigate weakly nonlinear hydrothermal two-dimensional convective flow in a horizontal aquifer layer with horizontal isothermal and rigid boundaries. We treat such a layer as a porous medium, where Darcy’s law holds, subjected to the conditions that the porous layer’s permeability and the thermal conductivity are variable in the vertical direction. This analysis is restricted to the case that the subsequent hydraulic resistivity and diffusivity have a small rate of change with respect to the vertical variable. Applying the weakly nonlinear approach, we derive various order systems and express their solutions. The solutions for convective flow quantities such as vertical velocity and the temperature that arise as the Rayleigh number exceeds its critical value are computed and presented in graphical form.


Introduction
A porous medium is a material that consists of a solid matrix with an interconnected void, and it is characterized by its porosity.Mathematical treatment and research on convection in porous media have been presented in many research articles.Darcy's law for porous media is the momentum equation, which is analogous to the Navier-Stokes equation.Heat transfer through a porous medium is a very common phenomenon.The natural tendency of fluid to expand when heated causes a density inversion to occur, if the heating is strong enough and a circulatory motion follows, termed convection.Convection in the fluid layer is a well-studied phenomenon, and it occurs in many natural settings: in the atmosphere, in the Earth's mantle and in many industrial applications including solidification and heating.Convective flows in a horizontal dendritic layer during alloy solidification are known to produce chimneys in the final form of the alloy.
The problem of thermal convection in an aquifer layer that occurs in many application areas, notably in environmental sciences and, in particular, in ground water flow [1,2], is fundamentally a problem within the category of convective flow in a porous medium.The effect of the aquifer nonhomogeneity on the onset of thermal convection, as well as on heat transfer through the aquifer has been analyzed [3].Many cases of convection through aquifers involve horizontal layers, which are non-homogeneous with variable permeability and thermal conductivity.
Riahi [4] studied convection when the boundaries have finite thermal conductivities.Vafai [5] investigated wall effects due to variable porosity.Riahi [6] used a perturbation approach for three-dimensional convection involving nonuniform temperature.He considered the case of a continuous finite bandwidth of modes [7].An aquifer model was presented by Fowler [8].Many studies have been presented by various authors [9][10][11][12][13][14][15] to investigate convection in porous media.The effect of variable permeability in porous media was analyzed analytically and numerically by Rees and Pop [16].Convection in cavities with conducting boundaries was carried out by Rees and Tyvand [17].The eigenvalue problem for cavities of various shapes was solved by them numerically.Rees and Barletta [18] described a linear stability analysis of a sloping porous layer with isoflux boundaries.Rees et al. [19,20] studied the effect of x-dependent variation in permeability.Bhatta et al. [21][22][23] studied weakly nonlinear convective flows in mushy layers with and without the magnetic effect and permeable mush-liquid interface.Rees and Barletta [24] analyzed the three-dimensional stability of a finite-amplitude convection in a porous layer heated from below.
The present paper studies the problem of nonlinear thermal convection in a horizontal aquifer layer with horizontal rigid and isothermal boundaries, and the aquifer layer is assumed to have variable permeability and thermal conductivity.The analysis is based on the weakly nonlinear theory [4] for the particular case of convection in a porous medium with rigid and isothermal boundaries and is extended here to porous layers with variable permeability and thermal diffusivity.We restrict our analysis to the non-dimensional form of the governing equations in a porous layer subjected to Darcy's law, and such a non-dimensional form of the governing system introduces three non-dimensional parameters, which are the Rayleigh number, hydraulic resistivity, which is the ratio of the viscosity to permeability, and non-dimensional thermal diffusivity, which is the ratio of diffusivity to a reference diffusivity at the lower boundary [3].The last two parameters are assumed to vary in the vertical direction, and the magnitude of the vertical rate of change of each of these two parameters is to be kept small.
We have found some interesting results, and in particular, we found that variations in the hydraulic resistivity or thermal diffusivity introduces asymmetry with respect to the mid-plane in the vertical profiles for the convective flow quantities.For the case where the vertical rate of change of either hydraulic resistivity or thermal diffusivity increases in the upward direction, then the flow appears to be stabilizing, while the flow is destabilizing if such a vertical rate of change decreases in the upward direction.

Governing Systems and Solutions
The non-dimensional system can be expressed as: Here, − → U , T , P,R, k, t, z, ξ, γ and ζ respectively represent the velocity, temperature, pressure, Rayleigh number, unit vector in the upward vertical direction, time, vertical coordinate, hydraulic resistivity, coefficient and diffusivity parameter.The coefficient γ is given by Rubin [3].
where ψ is the porosity, ρ is the fluid density, C s , C w are the specific heat of the solid and fluid, respectively, and ρ s denotes the density of the solid fraction.
For small δ (δ 1), we write: Here, ξ 1 and ζ 1 are constants.The boundary conditions are: where W is the vertical component of − → U .

Zeroth Order
To zeroth order in δ (to order 1 δ 0 , steady case), we have: Solutions are given by: where: and: Furthermore, we have:

Order of δ
The system in this case is given by: with boundary conditions: Multiplying (30) by φ 10 and (31) by − θ 10 R 00 , then we add them and average over the whole layer; and finally, applying (32), we obtain: Here, we define < .. > as the total horizontal and vertical average of the aquifer layer.Since: which yields: Using ( 29) and (34), we find the critical Rayleigh number R c at the onset of motion to be: It is seen from R c that its value is reduced for negative values of ξ 1 and ζ 1 , while the value of R c increases for positive values of these two parameters.Now, from (30) and (31), we can obtain: Let: Then, we have: To find the particular solution to (39) and (40), we take d 2 dz 2 − π 2 of (39) and use (40) to eliminate θ11 , which gives us: Assuming: We have: dz 2 = a 1 2 sin πz + 4πz cos πz − π 2 z 2 sin πz +a 2 2π cos πz − π 2 z sin πz − a 3 2π sin πz + π 2 z cos πz (44) and: Using the above in (42) and by setting both sides equal and setting similar coefficients to be equal, we obtain the following: Coefficients of z 2 sin πz yield: the coefficients of z sin πz, we obtain: The coefficients of z cos πz yield: The coefficients of sin πz yield: The coefficients of cos πz yield: Furthermore, we have: We use ( 43) and ( 45)-( 47) in (39) to find the particular solution for θ11 .where: Now, we obtain a complementary solution.Designating we can write: and using matrix algebra, (51) can be written in matrix form as: The solution to this homogeneous system is given by: Here, λ and K can be obtained from: det ( A − λI ) = 0 and ( A − λI ) K = 0 respectively.Now, det ( A − λI ) = 0 yields: Thus, we have four eigenvalues as: Let T be the eigenvector corresponding to eigenvalue λ i , i = 1, ..., 4. We have: Solving this system, we obtain: and we find four solutions: where the arbitrary constants C i i = 1, ..., 4 can be determined by using the four boundary conditions θ11 = φ11 = 0 at z = ± 1 2 .

Order of 2 δ 0
In this case, we have: We consider special solutions φ 10n , θ 10n of the linear system at order 1 δ 0 : where: Multiplying (57) by φ 10n and (58) by −R −1 00 θ 10n , adding them and averaging over the whole layer and, finally, applying (32), we obtain: Multiplying (62) by C n and taking the summation ∑ N n=−N , we find: Equations ( 57)-(59) give the solutions for φ 20 and θ 20 : where: For the two-dimensional case, (64)-( 67) become (since In this case, we have: Multiplying (70) by φ 10n and (71) by −R −1 00 θ 10n , adding them and averaging over the whole layer, we obtain: where: For two-dimensional flow, R 11 = 0.The present work presents the two-dimensional case, but we have provided the formulation for a general three-dimensional study so that it can be useful for future extension.
Next, the solvability condition at order 3 yields: Furthermore, the heat flux (heat transported by convection) can be estimated as follows: For the two-dimensional case: Using (34) in (77), we find that for the given value of R, the heat flux is higher for negative values of ξ 1 and ζ 1 , while the heat flux is lower for positive values of these parameters.

One-Dimensional Results for the Vertical Dependence of Linear Solutions:
In all our computations, the values of δ and used are 0.8 and 0.1, respectively.We have chosen 0.8 as the value of δ, even though it is assumed mathematically that δ 1, in order to make the effects due to variable hydraulic resistivity and diffusivity more visible in the graphs of the figures.Similar types of graphical considerations are made in other research areas.For example, Anderson and Worster [25] studied thin mushy layers with thickness δ 1, but in their graphical representations, they chose δ = 0.5 in order to show more visibly the effect of δ in their results.
The zeroth order solutions W 10 (z) and θ 10 (z) are displayed below in Figure 1. -  It can be seen from Figures 2 and 3 that vertical profiles for flow quantities to first order in δ, which are dependent on ξ 1 and ζ 1 , have lost their symmetry with respect to mid-plane flow and are stabilizing for ξ 1 > 0. The following graph in Figure 4   A similar asymmetry property holds here for Figures 4 and 5. Figures 6 and 7 present the linear solutions (W 10 (z) + δW 11 (z)) and (θ 10 (z) + δθ 11 (z)) with ξ 1 = 0.It can be seen from the presented linear results for the vertical dependence of the vertical velocity that for positive values of the vertical rate of change of the aquifer parameters (hydraulic resistivity and diffusivity), the magnitude of the vertical velocity decreases, which can imply flow stability, while an instability of the flow can result if such a vertical rate of change of the aquifer parameters is negative.4.2.Two-Dimensional Results for W 10 (x.z), θ 10 (x, z), W 11 (x, z), θ 11 (x, z), W 20 (x, z) : The two-dimensional zeroth order solution for W 10 (x, z) is shown in Figure 10, whereas Figure 11 displays the result for zeroth order solutions θ 10 (x, z).Two-dimensional views of flow quantities to zeroth order in δ show periodic behavior with period 2in the horizontal direction for a given z, while symmetric behavior in the vertical for a given x.
The results presented by these two figures show the same periodicity as in Figures 10 and 11, but with asymmetric behavior in the vertical for a given x.Similar behavior is evident in the next two Figures 14 and 15.
The two-dimensional profile of the nonlinear vertical velocity shows the same periodicity in the x-direction, and the vertical velocity shows slightly asymmetric behavior in the vertical direction for a given x value.

Conclusions
We investigated the problem of weakly nonlinear two-dimensional convective flow in a horizontal aquifer layer with horizontal isothermal and rigid boundaries.As in the application of such problems in the case of groundwater flow, we treated such a layer as a porous layer, where Darcy's law holds, subject to the conditions that that the porous layer's permeability and the thermal conductivity are variable in the vertical direction.The vertical axis is chosen positive in the upward direction anti-parallel to the force of gravity.In addition, our present study is restricted to the case that the subsequent hydraulic resistivity and diffusivity have small rates of change with respect to the vertical variable.We applied a weakly nonlinear approach, assuming a motionless and at most vertically-variable basic state and calculated the solutions for convective flow quantities such as vertical velocity and the temperature that arise as the Rayleigh number R exceeds its critical value R c .
We found that for the case where the vertical rate of either hydraulic resistivity or diffusivity is positive, then the convective flow is stabilizing, while the flow is destabilizing when the vertical rate of such quantities is negative.The results for the flow velocity indicated that asymmetry is projected due to variations of the hydraulic resistivity or thermal diffusivity.Our results also demonstrate the presence of a non-homogeneous effect in the aquifer layer, which is known to exist due to variations in permeability and/or thermal diffusivity.The heat transported by convection increases if the vertical rate of change of either hydraulic resistivity or thermal diffusivity decreases as z increases, while the heat flux due to convection decreases if the vertical rate of change of such quantities increases as z increases.
It should be noted that the present investigation of the aquifer layer was restricted, in particular, to the two-dimensional flow case.However, the present aquifer layer problem with small variations in thermal conductivity or permeability indicated possible flow structure in the form of hexagons, which will be left as one of our future subjects of investigations.In the case of multilayer porous media, we can explore in the future the possibility of the existence of stable hexagonal cells.