Convective Transport through Porous Layers

Convective motions are a multi-physics phenomenon, in which flow and transport processes interact in a two-way coupling. The density of the fluid depends on the value of transport variable and this back-coupling leads to non-linear behaviour. For the classical constellation of a closed fluid container heated from below convective motions appear, when a critical threshold for the Rayleigh number is exceeded. The heat transfer due to convection is much higher than in the case of pure conduction. Here systems of three layers are examined in detail. Using numerical CFD modelling it is shown that in layered systems different convective flow patterns appear than in the single layer case. The number and constellation of convection cells characterize steady flow patterns. Using a parametric sweep over the relevant parameter range of layer Rayleigh numbers and layer thicknesses we determine diagrams that show the excess heat or mass transfer of the dominant convection patterns, measured by the Nusselt- or Sherwood numbers.


INTRODUCTION
Heat and mass transfer through a system of porous layers are relevant topics in various constellations within our physical environment and in technical applications. Heat transfer in the sub-surface of the earth is very much determined by the involved geological layers. Aside from the understanding of natural geology, the heat flux is relevant in geotechnical systems as in heat storage beds [1]. The mass transfer of CO 2 is of concern for carbon sequestration in geological formations that are usually layered. Technical devices as packed-bed catalytic reactors contain layers [2]. Materials may consist of layers. When they are designed for thermal insulation, heat transfer is to be minimized [3]. Porous layers can be useful for maximizing heat dissipation, for example in electronic devices. In all these systems convection plays an important role.
Convective motions result from the coupling of flow and transport processes. They are thus a genuine multiphysics phenomenon. The nonlinearity of the coupling results in a remarkable behaviour of the studied systems. Here we deal with the classical constellation that a fluid of higher density overlies a fluid of lower density, which is discussed also as Horton-Rogers-Lapwood problem [4]. Above a critical threshold for the involved parameters, for example the density difference, a deviation from the trivial no-flow solution appears: steady flow patterns consisting of convection cells with circulating fluid can be observed. Following a classical analytical analysis, first rigorously presented by Lord Rayleigh for free fluids [5], the mutiphysics systems can be characterized by a dimensionless parameter combination, which nowadays is named the Rayleigh-number Ra or Ra-no. Below a critical value of Ra the system remains stagnant -above the critical value convective motions appear.
The mentioned coupling of flow and transport that leads to convection arises as effect of buoyancy caused by variable fluid density as the major coupling parameter. Fluid density mainly depends on temperature and salinity. Corresponding to this one speaks of thermal or haline convection. More complex behavior results in case of thermohaline convection when both temperature and salinity come into play. That additional complexity is not considered here.
The heat and mass transfer in a fluid system is significantly determined by the convective regime. It is convenient to describe the increased heat and mass transfer through the fluid layer as function of Ra. One may also characterize the transfer by dimensionless numbers: for heat transfer it is the Nusselt number Nu, for salt mass transfer the Sherwood number Sh. Both Nu and Sh are normalized to be 1 for the non-convective state with purely conductive heat or diffusive mass transfer.
In the classical description of convective motions, a single fluid layer is considered. Here the investigation is extended to systems of layers, through which heat or mass transfer appears due to diffusion and convection. The here-examined constellations consist of three horizontally aligned porous layers. However, they deliver clues for multi-layered systems in general. The general set-up, notation of geometrical entities and boundary conditions are depicted in Figure 1.
In order to restrict the parameter space in the three-layer system only constellations are considered in which the outer layers identical concerning the physical and geometrical properties, but differing from the intermediate layer. A dimensionless description is utilized, in which the geometry and the Ra-no.s are the only relevant parameters. As the porous medium Ra-no. is mainly determined by the permeability of the porous medium (see below), we may discuss the constellation in terms of aquifers and aquitards. In most cases the upper-and lowermost layers can be considered as aquifers, with a lower permeable layer sandwiched in between. If the intermediate layer is completely impermeable, convection rolls develop in the permeable layers above and below. If the sandwiched layer is only slightly less permeable (leaky aquifer) than the adjacent layers, convection rolls may develop that extend across the entire vertical extension of the system. By numerical modelling Bjørlykke et al. [6] examined the onset of convection in the very same set-up. The set-up, in which a permeable layer is sandwiched in between less permeable layers has attracted scientific interest, see Lipsey et al. [7] and references therein. Convective motions may appear in the sandwiched layer with an effect on the geothermal gradient within the layer. Here the study is extended with focus on a detailed examination of the flow patterns and the heat or mass transfer through such systems, covering the entire relevant parameter space.
The studied constellation differs from the partitioned porous layers, studied by Genç & Rees [8] and Rees [9], as here fluid flow in the aquitard is taken into account. It is a typical feature of geological layers that they are differently conductive to fluid flow as well as heat transport. The hydraulic conductivity of an aquitard can be several orders of magnitude lower than that of an aquifer. However, for moderate contrasts between the layers and long time scales the fluid flow within the leaky layer may not be ignored.
McKibbin & Tyvand [10] as well as Hewitt et al. [11] are dealing with systems in which thin very low permeable layers cut through permeable formations. In the limit case for very low permeabilities high Ra-no.s are required to induce convective motions. The systems that are in question here have comparably moderate dimensions and are thus far away from the limit case examined in the former studies.
The constellation with a porous layer overlain by a free fluid layer has attracted a lot of research studies [12,13,14,15,16,17], as there are several practical applications of interest. The topic here is a situation with several porous layers, governed by Darcy's Law. In the free fluid layer this approach is not valid. Instead such a layer is described by the Navier-Stokes equations. Although there convection patterns can be observed as well, the details concerning the onset of convection, cell shape and transfer are very different from the porous medium case that is in the focus here.

DIFFERENTIAL EQUATIONS & BOUNDARY CONDITIONS
Models of convective flow patterns are based on a system of partial differential equations, which result from mass and energy conservation formulations. Dependent variables are pressure, velocity and temperature or salinity, the latter depending on the type of convection, either thermal or saline. Relevant parameters are: the maximum density difference Δρ due to the different temperatures or salinities, the dynamic viscosity of the fluid μ, thermal or solute diffusivity D of the fluid-solid system and the permeability of the porous medium k.
Depending on conditions the resulting formulation can be simplified. The Boussinesq-Oberbeck assumption states that density variation is negligible except from its influence on buoyancy [18]. Another common assumption is that density is a linear function of either temperature (in the thermal case) or salinity (in the haline case).
The equations for the primitive variables pressure, temperature or salinity, can be transformed to a non-dimensional formulation with coefficients given by dimensionless numbers, which are combinations of the original parameters. For convection in porous media the porous medium Rayleigh number Ra is the only crucial characterizing dimensionless number: The porous medium Rayleigh number is defined in analogy to the Rayleigh number for free fluids. Details of the transformation to the non-dimensional formulation, briefly presented below, were given by Holzbecher [18] and more recently discussed by Dillon et al. [19] and Chandran et al. [20].
In the non-dimensional formulation, it is convenient to replace the pressure variable by the stream function . Using 2D Cartesian coordinates is implicitly defined by, where is the velocity field. denotes the horizontal and the vertical space direction. The formulation (2) ensures that the flow field is divergence free, i.e. , a condition that results from the mass conservation equation and the Boussinesq-Oberbeck assumption. The differential equation for is obtained using the vorticity vector , which in 3D is defined by the cross-product: Combining equations (2) and (3) for 2D flow in a plane yields only one nonzero vorticity component ω, which fulfills the equation: The dimensionless formulation is obtained, using a variable transformation with space unit H and time unit H2/D. H denotes the total height of the layered system. The velocity unit is thus D/H. By this non-dimensionalization the transport equation for heat or mass is given by, see [18]: is the normalized transport variable, i.e. temperature for thermal convection and salinity for haline convection. The terms in equation (5) represent the processes of storage, advection and diffusion. An explicit expression is obtained for vorticity: The sign in the vorticity equation of equation (6) is negative for the thermal case and positive for the for the haline case. The viscosity is assumed to be a constant, as we consider the same fluid with moderate changes of the transport variable only. In the numerical experiments reported here Ra mainly reflects the change of the permeability.
In order to complete the model formulation, boundary conditions have to be formulated. For a closed system the stream function has a constant value on all boundaries; without loss of generality we chose = 0. The conditions for at the outer boundaries are as follows: At the interfaces between the layers continuity conditions are required. The temperature and the heat transfer, in the thermal case, and salinity and mass flux, in the haline case, across the interfaces have to be identical in both layers that are separated by the interface. For both Ψ and θ this is achieved by setting a Dirichlet condition in one layer, and a Neumann condition in the other. If the variables in the adjacent layers are denoted by subscripts 1 and 2, the requirements for the streamfunction are: In the presented model for the upper-and lowermost layers we require the Dirichlet condition at the interfaces with the intermediate layer. For the intermediate layer the Neumann condition is demanded at both interfaces.
A similar construction for θ guarantees that diffusive heat and mass fluxes on both sides of the interfaces are identical. For the upper-and lowermost layers we require Neumann conditions at the concerned interfaces, and for the intermediate layer the demand is the Dirichlet condition.
In the following description we will concentrate on the thermal system, and not mention the haline case in particular. However, with the mentioned modifications concerning the sign in equation (6) and the outer boundary conditions all results are equally valid for the haline convection case as well.

SIMULATION
Numerical modeling is performed based on the system of equations (4)-(6). The coupled differential equations (4) and (5) are solved simultaneously. is computed using equation (2). ω is calculated using the explicit formula given in equation (6). The transient behaviour of the system is simulated. However, here we only discuss the steady states that are obtained after a sufficiently long simulation time.
The discretization is performed by Finite Elements using the software COMSOL Multiphysics [21]. COMSOL Multiphysics has been used for simulating convection patterns in several studies on porous media. Holzbecher [22] examined porous systems with open top boundary. Holzbecher [23,24] and recently Eckel & Pini [25] simulated flow patterns for very high Rayleigh-numbers in the context of CO2 storage. Systems with one horizontal interface have been investigated [26]. That investigation is extended here for the case of two interfaces.
Using COMSOL Multiphysics the differential equations (4) and (5) are treated using coefficient pde-modes. For the system with three layers altogether a system of 6 coupled differential equations has to be solved. For the discretization we choose finite elements on triangular or quadrilateral meshes with quadratic shape functions for all variables. The triangular meshes for the four geometries consist of 3516 and 3610 elements, which corresponds to approximately 15000 degrees of freedom (dof). In the parametric sweep we used quadrilateral meshes between between 6000 and 8000 elements (depending on d), which corresponds with up to 65526 dofs. In order to check mesh dependencies control runs are performed on even finer quadrilateral meshes with more than 22000 elements, corresponding to more than 180000 dofs.
For the solution of the transient simulation a time-stepping approach is used, with automatic timestep adjustment. The nonlinear equations are gathered in one system matrix. The resulting linear systems are solved by a direct solver.
As post-processing the Nusselt number Nu, for thermal convection (or Sherwood number Sh, in case of haline convection), are calculated. They represent the total heat (or mass) transfer through the system. in the non-dimensional formulation these numbers are defined by where the integral extends over an entire horizontal boundary. Nu (or Sh) respectively is computed using 4th order integration along the top or bottom boundary. In the steady state evaluations of equation (8) at the lower and upper horizontal boundary deliver the same value. This study focuses on steady state solutions. As the straight approach using the steady state versions of the equation (5) usually does not converge, we use long-term simulations of the unsteady equations to obtain the steady states! Only few of the simulations described below did not lead to a steady state. In that cases an oscillating convection pattern was reached instead.
The resulting flow patterns are generally not unique. They crucially depend on the initial condition. In the transient simulations for this study the initial values differ at a few spots from the linear pure diffusion profile. Without any such disturbance the computational simulation may not start, because the regime with no flow and linear profile is an analytical solution. At least one small disturbance is necessary to give a kick-start for the simulation, if it is physically unstable. The type of disturbance may have an influence on the final steady state. For example a single disturbance favours a system with two convection cells, a double disturbance favours a three cell pattern. However, despite of this the final state is not uniquely determined by the initial condition. A single initial disturbance may lead to steady states of 1, 4, 6 or 8 cells, as will be shown.
In some cases we observe an influence of numerical settings. As small local disturbances play a role in some of the investigated scenarios, it can be expected that changes of the numerical approach may have an effect on the simulation. We use triangular and quadrilateral meshes, which sometimes make a difference. In some cases different results are obtained, if the solution algorithm worked with scaled or non scaled variables.

Flow Patterns
The described model is run with a wide range of parameter variations. In order to limit the task, some constraints are set. It is assumed that the outer layers at the top and bottom of the model domain are of same type concerning their geophysical and geometrical properties. They have the same Ra-no. and the same thickness. Their Ra-no. in the sequel is referred to as Ra 1 , while the Ra-no. of the intermediate sandwiched layer is referred to as Ra 2 . In parametric sweeps Ra 1 varies between 100 and 500, Ra 2 varies between 10 and 300. The thickness of the intermediate layer d takes the dimensionless values of 0.1, 0.3, 0.5 and 0.7 from the total height H=1. Altogether these variations cover cases, where a low permeable aquitard is sandwiched in between thicker aquifers and cases, where a highly permeable layer lies in between less conductive layers.
For the different parameter combinations, some topologically distinct flow patterns can be observed. These are gathered and visualized in the following figures. Black horizontal lines depict layer interfaces. Moreover, the visualizations depict:  Decreasing the contrast of the layer characteristics, represented by the Ra-no.s, leads to cells with a higher curvature and the horizontal temperature gradients in the intermediate layer become more pronounced. Fig. 2b shows two convection cells, obtained for Ra 1 =125 in the aquifers and Ra 2 =90 in the sandwiched lower permeable formation. The latter type of pattern emerges quite frequently, especially in cases, in which the Ra-no.s take relatively nearby values.

Convection Modes
A parametric study is performed with the described model. The model is run for all combinations of parameters, listed in Table 1. The range of the parameters represents the relevant cases with steady convection flow patterns. For higher Ra-no.s one partially enters in the region, where oscillatory flow pattern become more likely. A transition to oscillatory patterns was observed in some model runs here. For single homogeneous layers Holzbecher [27] investigated this transition, but this is not in the focus here. For each of the four variations of the geometry we make a parametric sweep with variations of Ra 1 and Ra 2 , i.e. 132 simulations. In all cases a steady state flow is found, that falls into one of the pattern types described above. For few parameter combinations the cell pattern is not unique: depending on the initial conditions and options for the numerical solution different modes are obtained. In those cases several convection solutions exist. For few combinations also oscillating patterns are observed. For all runs in this study, for each parameter set at least one steady solution emerges. Figure 3 shows the dominant cell pattern observed for the four geometry constellations with different layer thicknesses. The dominant patterns are those with a single cell, as shown in Fig. 2e, with two cells, as shown in Fig.s 2a and 2b, and with four cells, as shown in Fig.s  2c and 2d. The 0-cell cases represent no convective, i.e. pure diffusive flow.
In case of the thin intermediate layer (d=0.1) the 4-cell mode is observed in almost all cases, where an aquitard is sandwiched in between two aquifers. If the intermediate layer has a lower Ra-no. the 2-cell mode prevails. For Ra1=100 the single cell mode can be preferred. The flow modes that do not appear in Figure 3, i.e. the ones depicted in Fig.s 2e-2j, are observed for relatively few cases of the parametric sweep. Depending on the chosen initial conditions and numerical options, as outlined above, for the same parameter combination also one of the dominant modes, considered in Fig. 3, appears. For d=0.3 we find Ra-No. combinations leading to the pattern, shown in Fig. 2g and to the 6-cell mode.

Heat-or Mass Transfer
The Nusselt-no. Nu, defined by equation (8), represents the heat transfer through the layered system. For the pure diffusive situation with no flow in the entire domain Nu is 1. Nusseltno.s above 1 measure the increased transfer due to convection.
Nu depends on the layer Ra-no.s and on the geometry. For the here considered model setup the geometry is uniquely represented by the thickness d of the intermediate layer. Moreover the transfer is also dependent on the convection mode. For the dominant modes of the parametric sweep the Nu-numbers are depicted in Figure 4. For all geometries the range of the Nu-no.s is between 1 and 7.
For the thin intermediate layer (d=0.1) Nu-no. results show a strong dependence on Ra1, which characterizes the higher permeable regions. Between 150<Ra1<300 the contour lines are almost vertical, indicating no dependence on Ra2. The transition to higher transfer rates for Ra1=100 and high values of Ra2 can be explained by a change of the convection pattern. As the corresponding sub-plot in Figure 3 shows, the 1-cell mode is dominant for that parameter combination. For d=0.3 a small irregularity of the prevailing pattern can be observed in the lower left corner and a strong deviation on the lower right corner of the figure. Both emerge from a change of the dominant flow pattern, as can be seen in the corresponding sub-plot in Fig. 3.
For d=0.5 there is a stronger dependency on Ra2 in the lower part of the plot, but a stronger dependency on Ra1 in the upper part. A comparison with Figure 3 shows that that roughly coincides with the modes. In the part with dominant 4-cell mode Ra2 is more important, while for the 2-cell pattern Ra1 plays a bigger role.
For a thick intermediate layer, d=0.7, the just described behavior is still recognizable, but on a weaker scale. Comparison with Figure 3 shows that the 2-cell pattern is the dominant mode for this geometry.

DISCUSSION & CONCLUSION
In a layered system convection patterns appear that cannot be observed in homogeneous systems. Convection cells emerge that are connected to layers: there are solutions with 2x2, 2x3 and even 2x4 cells in the square unit domain. Moreover superposition of cells occurs, when within a larger cell two smaller cells appear. Various complex circulation patterns emerge that are described above.
Heat and mass transfer depends on the geometrical constellation, i.e. here the height of the different layers and the Ra-numbers, as given parameters. It also depends on the flow pattern, which is not unique for the non-linear system. This can be of interest in technical systems, in which minimal or maximal transfer is to be achieved. The current study gives an aid on which parameters can be chosen do reach a state of optimum heat or mass transfer.
Attention is to be paid on the decrease of heat transfer with Ra 2 in the high Ra-numbers regions. It can be observed in the upper left corners of sub-plots in Fig. 4 for d=0.1, 0.3 and 0.5. Nu can generally be expected to increase with increasing Ra, what the mentioned figures seem to contradict. However, a comparison with the corresponding sub-plots in Fig. 3 reveals, that there is a transition of the flow-pattern in the concerned regions of the plots. Raising Ra 2 there can lead to a change from the dominant 2-cell pattern to the 4-cell pattern. This is accompanied with a change of the Nu-no., which is lower for the 4-cell mode. However, detailed investigation of the obtained results shows that the mentioned decrease of Nu-no.s appears also without mode change -a phenomenon that is not completely understood.