Density-driven convection in an inhomogeneous geothermal reservoir

International Journal of Heat and Mass Transfer 127 (2018) 784–798

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Density-driven convection in an inhomogeneous geothermal reservoir E.B. Soboleva Ishlinsky Institute for Problems in Mechanics, Russian Academy of Sciences, 119526 Moscow, pr. Vernadskogo 101, b. 1, Russia

a r t i c l e

i n f o

Article history: Received 10 May 2018 Received in revised form 12 July 2018 Accepted 6 August 2018

Keywords: Geothermal system Porous domain Admixture transport Porosity Permeability Salinity Diffusivity Darcy’s law Haline convection Numerical simulation

a b s t r a c t Present study is motivated by problems of admixture transport in geothermal systems. Numerical simulations of stochastic haline convection in an inhomogeneous domain consisting of a homogeneous isotropic porous medium divided by a horizontal layer at different porosity and permeability are performed. Convection is driven by more dense brine near the upper boundary which is fed by a dissolved admixture diffusing from the boundary. An influence of interior layer on convective flows and mass transfer in all domain is investigated. As obtained, dynamic behavior of fluid in the layer is forced or natural mode of convection depending on properties of this layer. There is the intermediate case when forced convection arisen initially is transformed into natural convection with time. In this case, the layer can order convective flows above itself so that the roughly periodic structure of ‘‘brine balls” is formed. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction The problem of admixture transport in groundwater saturating soils and rocks has a fundamental importance for alternative power engineering, ecology, environmental management and other human activities. In many cases, transport of dissolved admixture in the Earth is caused by the temperature and salinity gradients leading, in the gravity field, to a growth of instability and development of natural convection [1,2]. In deep aquifers, salinity gradients can be formed by a salt deposit nearby or by brine flows carrying dissolved admixture. There are several geological reasons inducing temperature inhomogeneities such as magmatic heat sources, chemical reactions or geothermal gradient which is resulted from the temperature increase within the Earth in the direction to the hot Core and available all over the place [3]. Normally, the geothermal gradient is about 20–30 K per each kilometer. Estimations exhibit that, in real geological conditions, convective motions of concentrated brine due to the geothermal gradient may be appreciably less than due to the salinity gradient [4]. In these conditions, thermal convection can be disregarded and only haline convection be taken into consideration. Haline or, as is called also, density-driven convection is widely discussed in the context of carbon capture and sequestration techE-mail address: [email protected] https://doi.org/10.1016/j.ijheatmasstransfer.2018.08.019 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved.

nology, wherein liquid carbon dioxide is injected into geological formations to store it underground [5]. Diffusion of carbon dioxide and its dissolution in saline aquifers leads to a local increase in the brine dense and development of downward convective motions (see [6–13] and references herein). Similar problems on convection are actual in respect to exploiting underground geothermal systems which use deep heat of Earth and are alternative renewable energy sources [14,15]. If there is a salt supply near a geothermal reservoir, haline convection can develop and enhance an admixture transport into hot water to be extracted. There are quite a lot of publications on numerical study of density-driven convection in rocks arisen from a horizontal unstable denser layer of brain. In these works, they investigate nonlinear convective regimes in a comparison with a linear stability analysis [6], an influence of the Rayleigh number on the intensity of transient convection in an initial stage [7], effects of random fluctuations of porosity and permeability on the onset of convection [8], mass transfer quantities based on model experiments and on comparable numerical simulation [10], convective regimes from onset to shutdown in anisotropic media having different values of horizontal and vertical permeability [11], an interplay with forced [16] or thermal modes of convection [12], and a complicated case including geochemical reactions [13]. One of key factors defining parameters of convective flows and mass transfer is the structure of solid matrix containing a liquid

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

785

Nomenclature aspect ratio, ¼ L=H [–] concentration (mass fraction) of admixture, ¼ qc =q [–] solubility (mass fraction of admixture in a saturated sat solution), ¼ qsat [–] c =q D coefficient of diffusion in a porous medium [m2 s1] Dc coefficient of molecular diffusion for a free solution [m2 s1] e unit vector in the direction of gravity acceleration, =(0, 1) [–] g module of gravity acceleration [m s2] H height of porous domain [m] hd distance from the upper boundary of domain to interior porous layer [m] hl width of interior porous layer [m] depth of admixture penetration [m] hp k permeability [m2] L length of porous domain [m] M dimensionless total admixture mass per the unit horizontal length [–] M0 dimensionless total admixture mass per the unit horizontal length due to diffusive mass transfer [–] n unit vector normal to boundary, [–] P pressure [Pa] Rm relative admixture mass, ¼ M=M 0 [–] Rd Rayleigh-Darcy number [–] S dimensionless density [–] t time, except Sections 2.3 and 4.3 [s], [d]; dimensionless time, in Sections 2.3 and 4.3 [–] Dt numerical time step [–] U ¼ ðU x ; U y Þ vector of filtration velocity ¼ /W [m s1] u ¼ ðux ; uy Þ vector of dimensionless filtration velocity, ¼ /w [–] W ¼ ðW x ; W y Þ vector of fluid velocity [m s1] w ¼ ðwx ; wy Þ vector of dimensionless fluid velocity [–] jWjmax maximal module of fluid velocity [m s1], [mm d1] jwjmax maximal module of dimensionless fluid velocity [–]

A c csat

phase. In certain Earth’s zones, a layered structure of crust has been formed under the high pressure during different geological events therefore a consideration of horizontal layered porous media is the way to approach real conditions. Investigators consider systems consisting of several separately homogeneous isotropic layers or anisotropic systems with continuous variations in porosity and permeability as an approach to multilayered porous media [17– 20]. However, the data showing how a single interior porous layer at different porosity and permeability influences haline convection are not enough. There are only little results on this problem based on a simplified mathematical model and presenting some integral characteristics of initial stage of convection [4]. In the present paper, we continue investigations starting with [4]. The objective of our work is to study stochastic haline convection from the beginning to a developed stage in a porous domain containing a horizontal interior porous layer with lower porosity and permeability. Convection is driven by an unstable dense fluid layer near the top of domain, which is formed due to admixture diffusion from the upper boundary of this domain. The paper is organized as follows. In Section 2, we describe the problem under study and discuss the mathematical model suitable for solving this problem. In contrast to [4], we discuss how the permeability and coefficient of diffusion relate to the porosity and introduce appropriate relations into our description. The physical parameters of considered geothermal system based on reference data are put in Section 3. A brief report on the finite-difference numerical code designed for solving 2D problems is given in Section 4. We write

x

Dx y Dy1 Dym

horizontal coordinate, except Section 2.3 [m]; dimensionless horizontal coordinate, in Section 2.3 [–] numerical space step in the x-direction [–] vertical coordinate, except Sections 2.3 and 4.3 [m]; dimensionless vertical coordinate, in Sections 2.3 and 4.3 [–] maximal numerical space step in the y-direction [–] minimal numerical space step in the y-direction [–]

Greek symbols salinity expansion coefficient, =0.815 [–] dimensionless permeability [–] dynamic viscosity P dimensionless pressure [–] q density of solution [kg m3] qsat density of saturated solution [kg m3] qc density of dissolved admixture [kg m3] qsat density of dissolved admixture in a saturated solution c [kg m3] q0 density of pure water [kg m3] r dimensionless coefficient of diffusion [–] s tortuosity factor [–] / porosity [–]

a j ls

Superscripts ⁄ upper boundary of domain Subscripts 1 main porous medium, except Dy1 in Section 4.3 2 interior porous layer x horizontal axis y vertical axis th threshold value

also about tests and validation of the code. Local and integral numerical quantities characterizing flows and mass transfer are defined in this section as well. In Section 5, the report on numerical simulations of haline convection is given. We starts with a basic case corresponding to convection in a homogeneous domain with no interior layer. Next, the problem in a domain containing a horizontal layer at different properties is solved. We study how the layer properties influence selecting the mode of convection in the layer and convection characteristics in a domain as a whole. Effect of location of layer relative to the upper boundary and its width are discussed as well. We summarize our results and draw a conclusion in Section 6.

2. Theoretical background 2.1. Problem under study A rectangular domain being a part of reservoir infinite in the horizontal direction is considered. The domain consists of homogeneous isotropic porous medium with uniform porosity and permeability (referred to a main medium) everywhere except an interior horizontal layer which is at different porosity and permeability. The interior layer is located at the distance hd from the top and has the width hl . Initially, the reservoir is filled with pure motionless water. The constant concentration of admixture is held at the upper boundary modeling an admixture source. No-flux conditions

786

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

are applied to the vertical boundaries of domain. The fluid velocity is equal to zero at the horizontal boundaries of domain. The sketch of the problem is given in Fig. 1. Initially, an admixture from the top diffuses into water forming a brine layer near the top. As the density of solute increases with admixture concentration, denser fluid near the top is found to overlie less dense fluid. When the admixture concentration in water reaches a certain value, the profile of solute density becomes unstable in the gravity field and try to fall. Small disturbances of velocity grow leading to a development of natural density-driven convection. Note, the instability of fluid can develop independently of gravity, for example, under the capillary pressure gradient [21,22]. However, such effects are not taken into consideration in the present study. The objective of our work is to investigate how the interior layer at different properties affects convective motions and mass transfer in the reservoir. We stop our study before convective flows reach the bottom of reservoir so that the problem defined approaches the problem in a semi-infinite space. 2.2. Mathematical model We use the continuity, Darcy’s and admixture transport equations. Water is assumed to be incompressible. The density of solute varies only with the mass of dissolved admixture and is governed by the linear equation of state. The governing equations are as follows [1,2]

rU¼0 U¼

/

ð1Þ

j ðrP  qg  eÞ ls

@ qc þ U  rqc ¼ rð/Drqc Þ @t

q ¼ q0 þ aqc

ð2Þ ð3Þ ð4Þ

The pressure is subjected to the hydrostatic effect. The pressure distribution at the initial instant Pin follows from the Darcy’s equation at U ¼ 0. The density is assumed to be constant initially: q ¼ q0 . We have

rPin ¼ q0 g  e

ð5Þ

Integrating Eq. (5), a linear equation showing a decrease in the pressure with the height is obtained. The initial conditions have the form

t ¼ 0 : U ¼ 0; q0 ¼ 0; P ¼ Pin ¼ P  þ ðP b  P ÞðH  yÞ=H

ð6Þ

Here, P  and Pb are the pressure at the top and bottom of considered domain. The boundary conditions are

x ¼ 0; x ¼ L : n  U ¼ 0; n  rP ¼ 0; n  rqc ¼ 0

y ¼ H : n  U ¼ 0; P ¼ Pb ; n  rqc ¼ 0 y ¼ 0 : n  U ¼ 0; P ¼ P ; qc ¼ qc

ð7Þ

One should remark on the model employed. The permeability of porous medium k in the Darcy’s Eq. (2) correlates with its porosity, pore size, structure and other physical properties. The classical relation is the Kozeny-Carman equation used extensively. This equation can be written for k as a function of porosity / with the coefficient C k to include an influence of other factors [1]

k ¼ Ck

/3 ð1  /Þ2

ð8Þ

The coefficient C k is determined explicitly and the relation (8) is modified in special cases [23–26]. However, in our consideration, we do not specify intrinsic properties of porous medium and accept Eq. (8) with a constant C k as the basic relation. The right-hand side of the equation of admixture transport (3) includes, in a common case, the hydrodynamic dispersion D to be a second rank tensor [2] instead of the scalar D. The tensor D is introduced taking into account a finiteness and tortuosity of porous space and can be divided into two parts, one of these is the mechanical dispersion Dm and the other is the effective diffusion Dc T. A second rank symmetric tensor T represents the tortuosity of porous medium whereas Dc is the molecular diffusion for a free solution. Determination of Dm and its effect on transport phenomena are discussed in many works [2,27–32]. In our paper, we consider intensive convective flows with the advective flux to be major in a comparison with the dispersive flux and therefore, in the present level of study, neglect the mechanical dispersion. However, it is important to take into account the effective diffusion because an admixture from the upper boundary spreads down into water due to diffusion mechanism. In an isotropic porous medium, the components of T are represented as sdi;j where s is a scalar referred to the tortuosity factor and di;j the Kronecker delta. Under accepted assumptions, the dispersion coefficient in the admixture transport Eq. (3) reduces to the form

D ¼ sD c

ð9Þ

The tortuosity factor s is determined in different empirical models which, in saturated porous media, define s depending on the porosity / [24,26,33,34]. For example, the first proposed relation s ¼ /4=3 [35] is suitable for rough modelling in many cases. However, this relation was derived for gaseous diffusion in soils and, at small /, gives appreciably lower values of s as it follows from experiments. We consider porous media at low porosity / ¼ 0:05  0:2 and employ the relation s ¼ /0:4 which provides a good agreement with experimental data at / < 0:5 [33] and agrees with theoretical predictions [26]. So, we have

D ¼ /0:4 Dc

ð10Þ

The relations (8), (10) allow one to consider the permeability and effective coefficient of diffusion as a function of only porosity. The equation of state (4) includes the salinity expansion coefficient a which is a constant. 2.3. Dimensionless form of the governing equations First, the Darcy’s equation is modified by introducing the deviations of density and pressure from their initial values q  q0 ¼ aqc and P  Pin , respectively. Subtracting Eq. (5) from Eq. (2) and performing simple transformations, the Darcy’s equation reduces to the form

U¼ Fig. 1. Sketch of the problem.

j ðrðP  Pin Þ  aqc g  eÞ ls

ð11Þ

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

787

The governing Eqs. (1), (2), (11) in terms of U, P  P in , qc are taken into consideration. We put these equations into dimensionless form.

3. Physical parameters of geothermal reservoir

The height of domain H, diffusive velocity Dc =H, time H2 =Dc , density

The temperature of Earth’s crust increasing due to the geothermal gradient reaches 450–550 K in the depth of several kilometers. The pressure can reach tens-hundreds atmospheres. We consider a geothermal reservoir at T = 500 K and P = 1.5  107 Pa. At these conditions, the parameters of water and dissolved admixture are as follows

sat aqsat c , and pressure aqc gH are used as scales. According to the equation of state, the density scale aqsat c is equal to the density difference sat  q0 . The permeability and of saturated brine and water, aqsat c ¼ q

coefficient of diffusion are normalized by the reference values k1 , and Dc , respectively. The dimensionless velocity u ¼ U H=Dc , pressat sure P ¼ ðP  P in Þ=ðaqsat c gHÞ and density S ¼ qc =qc are considered. The set of dimensionless equations is the following

ru¼0

ð12Þ

u ¼ Rd jðrP  S  eÞ

ð13Þ

q0 ¼ 842:5kgm3 ; ls ¼ 1:212  104 Pas; csat ¼ 0:3417; Dc ¼ 3  109 m2s1 ; a ¼ 0:815

ð21Þ

The boundary conditions are transformed into equalities

The magnitudes of q0 and ls are defined in [36]. The data on csat and Dc are presented in [37]. The salinity expansion coefficient a follows from to the linear approximation of sea water equation of state [38]. The porosity and permeability of porous medium in geothermal conditions are really variable, so that the first is typically about 0.2 and less, and the second is from 1017 m2 to 1012 m2. The effective coefficient of diffusion is calculated by Eq. (10). For the main porous medium, we put

x ¼ 0; x ¼ A : n  u ¼ 0; n  rP ¼ 0; n  rS ¼ 0

/1 ¼ 0:2; k1 ¼ 1013 m2 ; D1 ¼ 1:576  109 m2 s1

y ¼ 1 : n  u ¼ 0; P ¼ 0; n  rS ¼ 0

At the upper boundary, we put the concentration of dissolved admixture c

/

@S þ u  rS ¼ rðr /rSÞ @t

ð14Þ

The initial conditions reduce to the form

t ¼ 0 : u ¼ 0; S ¼ 0; P ¼ 0

y ¼ 0 : n  u ¼ 0; P ¼ 0; S ¼ S

ð15Þ

ð16Þ

The parameters in Eqs. (12)-(14), (16) are the Rayleigh-Darcy number, Rd, the normalized permeability and coefficient of diffusion, and the density of admixture at the top

Rd ¼

aqsat c gk1 H ; j ¼ k=k1 ; r ¼ D=Dc ; S ls Dc

ð17Þ

The values j, r are assumed to be j1 , r1 in the main medium whereas j2 , r2 in the interior layer. The Rd number does not characterize the convective process completely because includes the maximal admixture density qsat c and coefficient of molecular diffusion Dc rather than the values driving convection really. Considering the main porous medium, convective flows depend on the density at the top qc and effective coefficient D1 therefore the  actual Rayleigh-Darcy number Rd is defined as follows 

Rd ¼

agk1 Hqc S ¼ Rd ls D 1 r1

ð18Þ



The Rd number is the parameter of convection similarity. To quantify the amount of admixture dissolved in water, we use the concentration c as the most common parameter interchangeably with the dimensionless density S. These values are related as

S¼

ð1  acsat Þc csat ð1  acÞ

ð19Þ

c ¼ 4:0  102

ð23Þ

The space sizes of considered domain are

H ¼ 25m; L ¼ 50m

ð24Þ

The properties of interior layer depending on the porosity /2 vary in computations with /2 . The position and width of interior layer determined by hd and hl vary as well. 4. Numerical code, tests and quantities 4.1. Method and code We employ the finite-difference 2D numerical code designed for solving the problems of haline convection in porous media. The earlier version of this code was successively used for numerical simulation of convective flows and mass transfer in groundwater resulted from the salt accumulation near the front of phase transition due to the upward solute flow [39,40]. Later, the code was used for some simulations of admixture transport in geothermal reservoirs during an initial stage [4]. The governing differential equations written in Cartesian coordinates are approximated by finite-difference equations on a straggled non-uniform grid (Fig. 2). Scalar variables are defined in the centres of meshes whereas components of velocity are at their

At the boundaries of interior layer, the pressure and density of admixture are continuous. The component of filtration velocity normal to the boundary is continuous as well [18]. It means that, in the case of horizontal layer, the vertical component of filtration velocity does not change whereas the horizontal component changes by step at interfaces. As a result, flow lines cross the boundary between a main medium and interior layer with changing in flow direction; this effect is discussed in Section 5.2. The fluid velocity is discontinuous because of changing in the porosity. The conditions at the boundaries of interior layer are

P2 ¼ P1 ; S2 ¼ S1 ; ux2 ¼ ux1 j2 ;uy2 ¼ uy1 ; wx2 ¼ wx1 j2

ð22Þ

/1 / ;wy2 ¼ wy1 1 /2 /2 ð20Þ Fig. 2. Space grid.

788

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

faces. The calculation procedure consists of several steps. First, we use the SIMPLE-type algorithm [41] to define the coupled velocity and pressure from the finite-difference analogues of Eqs. (12) and (13). Second, we integrate the admixture transport Eq. (14) written in a finite-difference form and obtain the density. Finally, corrections of all variables are performed by an iterative scheme. We use two-point and three point central-difference schemes for first and second partial derivatives, respectively. The discretization procedure gives rise to systems of simultaneous algebraic equations tridiagonal in x- and in y-directions. Each algebraic system is solved by Thomas’ algorithm applied successively in x- and ydirections [42]. Note, finite-element and spectral methods give more accurate results for convection problems comparing with the finitedifference method [43]. However, we choose the finitedifference discretization because of its ease of coding, flexibility and economy. A high accuracy of our solutions is achieved by employing a fine space grid and appropriate time step. In addition, our discretization procedure with the use of rectangular straggled space grid leads to approximate algebraic equations identical to those in the finite-volume method satisfying conservation principles. 4.2. Numerical quantities One can estimate dimensionless parameters of convection in a homogeneous domain from Eqs. (17)–(19) at physical properties set in Section 3. We have

Rd ¼ 2:194  104 ; j1 ¼ 1; r1 ¼ 0:525; 

S ¼ 8:731  102 ; Rd ¼ 3:647  103

ð25Þ 

The actual Rayleigh-Darcy number Rd is in two orders of magnitude higher than the threshold value Rdth = 4p2  40 corresponding to the onset of convection in the concentration RayleighBenard problem (or the Horton-Rogers-Lapwood problem at the  zeroth temperature gradient) [1]. As Rd is so high, we predict a development of stochastic haline convection which can be treated with the use of some local, average and summary quantities. We try to avoid an effect of vertical boundaries and approach computed quantities to ones in a domain infinite in the horizontal direction. For this reason, in relations (26)-(28), (30) below, we summarize in x-direction over internal meshes from ni to ne (see Fig. 2) excluding some near-boundary areas from our summarizing procedure. As typical, we put ni and ne so that to exclude the areas of length of 0.2 at the left and at the right. The length of domain taking into account for summarizing procedure is lx. We find the maximal fluid velocity jwjmax from the filtration velocity field as follows

0

jwjmax

1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi B1 C u2xi;j þ u2yi;j ji ¼ n ; n A ¼ max @ e i /i;j j ¼ 2; m

An admixture penetrates into a domain only through the upper boundary due to diffusion therefore the global mass balance assumes that the admixture mass M in the reservoir at any time te should be equal to the amount of admixture M J diffusing through the upper boundary from the beginning till this time t e . Let the mass flux through the boundary at the random time t be denoted as J k . The total mass balance is written as

M ¼ MJ ; MJ ¼

X

J k Dt; J k ¼ /1 r1

k¼1;ke

1 X Si;mþ1  Si;m Dxi lx i¼n ;n Dymþ1=2 i

ð28Þ

e

Calculating M J , the mass flux J k multiplied by the time step Dt is summarised over all time steps from 1 to ke corresponding to t e . However, numerical errors during computations lead to some discrepancy in the first Eq. (28). One can evaluate an accuracy of computations by examining the first Eq. (28). We estimate an error dM at each time level. From the definition,

  M  MJ   dM ¼  M þ MJ 

ð29Þ

A high level of accuracy of numerical solutions is verified by a little value of dM not exceeding 1% in all computations. The relative admixture mass Rm is considered to estimate a convection contribution into enhancing mass transfer. The value of Rm is calculated as Rm ¼ M=M 0 . Here, M is the admixture mass defined above whereas M 0 is the admixture mass in a model problem at g ¼ 0 that is at zero gravity conditions to exclude convective motions. We obtain M and M 0 sequentially and then calculate Rm . The mass transfer is characterised by the depth of admixture penetration Hp . We calculate initially the average admixture density qj per the unit horizontal length defined at the level of coordinate y or, in a discrete form, at the space grid number j corresponding to y.

qj ¼

1 X Si;j /i;j Dxi lx i¼n ;n i

ð30Þ

e

As obvious, qj is maximal at the upper boundary (y ¼ 0) and decreases with distancing from it. We define the depth of penetration hp as a distance from the upper boundary to the level of y at which the appropriate value qj is equal to 0.01 fraction of q ; the value of q is calculated at y ¼ 0. So,

qj =q ¼ 0:01

!

hp ¼ jyj

We multiply hp by the scale and present below in a dimensional form. Note, the admixture mass M can be obtained by summarizing P qj as M ¼ j¼2;m qj Dyj .

ð26Þ

The components of filtration velocity ux and uy are defined in the centre of (i,j) mesh by interpolation of values at mesh faces. For the porosity, the equality /i;j ¼ /2 is used if the grid node (i,j) locates in the interior layer, otherwise, /i;j ¼ /1 . The value of jwjmax is transformed to a dimensional form and given below as jWjmax in units of millimeter per day. The total admixture mass per the unit horizontal length M is found by summarizing the local density multiplied by porosity over space meshes Dxi and Dyj :

M¼

1X X Si;j /i;j Dxi Dyj lx j¼2;m i¼n ;n i

e

ð27Þ

Fig. 3. Diffusion problem: comparison of analytical (Eq. (34)) and numerical solutions.

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

4.3. Validation of the code For accurate computations of stochastic haline convection, we need to use fine space grids. After series of testing, we have chosen for major computations the 5001  2451 grid uniform in the

789

x-direction and non-uniform in the y-direction with reducing in the step size near the upper boundary. The step is Dxi = Dx = 4.0  104 in the x-direction. In the y-direction, it varies from Dy1 = 4.24  104 at the bottom to Dym = 4.24  105 at the top. The time step is typically Dt = 8.64  109 corresponding to 30 min of real time. Numerical study has been performed with the use of PC with 4-core processor Intel Core i7-920. Nearly 80 h of computation time or less are spent on each calculation. The space and time steps are satisfied by certain requirements. First, one should use a refined grid near the upper boundary to resolve accurately an admixture entry into a domain due to diffupffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi sion. An admixture diffuses over the distance sDt ¼ /0:4 Dt durqffiffiffiffiffiffiffiffiffiffiffiffiffiffi ing the time step Dt, therefore we require Dym < /0:4 Dt to have grid points in the diffusion boundary layer form the beginning of qffiffiffiffiffiffiffiffiffiffiffiffiffiffi computations. We obtain /0:4 Dt = 6.74  105 and really Dym < 6.74  105. The numerical code is examined by two problems solved on the same space grid and using the same time step as declared above. The test problems are in zero-gravity conditions (g ¼ 0) and reduced to 1D consideration. To validate a high level of accuracy of admixture diffusion simulation, we compare our numerical result with an analytical solution. At zeroth gravity, fluid is motionless and the admixture transport Eq. (14) reduces to the diffusion equation. One can find the analytical solution of the diffusion equation in the case of uniform porosity and coefficient of diffusion. If an admixture spreads down in a semi-infinite space from the boundary y ¼ 0 kept at the constant S , the analytical solution is as follows [44]

Fig. 4. Concentration step problem: comparison of moving step with a static step located in the level of y = 0.25 (a), 0.50 (b).

   y Sðy; tÞ ¼ S 1  erf pffiffiffiffiffiffi 2 rt

Fig. 5. Field of admixture concentration in the initial stage of convection at the time instants t = 200 (a), 245 (b), 312 (c) days.

ð31Þ

790

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

Here, erf ðzÞ ¼ p2ffiffipffi

Rz 0

2

en dn is the error function. The amount of

admixture M 0 ðtÞ per the unit boundary length spreading into the volume by the time t is found by integrating the mass flux r @S=@yjy¼0 multiplied by / over the time. We have

Z M 0 ðtÞ ¼ / r

t 0

@S=@yjy¼0 dH

ð32Þ

To integrate Eq. (32), one should substitute the function Sðy; tÞ   yffiffiffiffi is represented by the Taylor from Eq. (31). The function erf 2p rt pffiffiffiffi 3 series as erf ðzÞ ¼ 2= pðz  1=3z þ . . .Þ [45]. Because of the vari  yffiffiffiffi to be small at y ! 0, we hold only the first term of able 2p rt the series and reduce the equation on M 0 ðtÞ to the approximate equation which after simple transformations takes the form

1 M 0 ðtÞ ¼ / r S pffiffiffiffi

p

Z

t

0

1 pffiffiffiffiffiffiffiffi dH

rH

ð33Þ

Integrating Eq. (33), we find the explicit dependence M 0 ðtÞ as follows

M 0 ðtÞ ¼

pffiffiffiffi 2/ S r 1=2 pffiffiffiffi t

p

Fig. 6. Field of admixture concentration in the developed stage of convection at the time instants t = 500 (a), 1275 (b) days.

ð34Þ

The relation (34) written for the main medium at /1 , r1 and calculated at parameters (22), (25) is exhibited in Fig. 3. For a comparison, we solve numerically the problem in a homogeneous domain at g ¼ 0 and other parameters being the same. The numerical solution on M 0 ðtÞ is plotted in Fig. 3 as well. As clear, the analytical and numerical data are in excellent agreement showing a highly accurate simulation of admixture diffusion. We also test our numerical code by the problem of concentration step moving at uniform velocity. We examine a simulation of high concentration gradients occurring in intensive convective flows. In this problem, pure water moves downward at the velocity wy = 4  103. The upper boundary is hold at the constant admixture density S , so that, starting from the initial time, the brine at S enters a domain and forms a concentration step. This step moves down with a forced flow and spreads due to diffusion. We set g ¼ 0 to exclude natural convection and keep parallel moving. The velocity of forced flow wy is chosen to be close to the magnitudes of convective velocity in numerical solutions of the problem under study given below (see Fig. 7(a)). In Fig. 4, the fragment of concentration profile at two different times is shown; the zone of high concentration gradient is displayed. For a comparison, we place the motionlessness concentration step at the level y = 0.25 (a) or y = 0.50 (b) and solve numerically the problem at wy = 0. We obtain spreading of admixture due to only diffusion not disturbed by a flow. The marks in Fig. 4 are plotted in three grid point intervals. As is clear, the downward flow in the first case does not distort the density profile: deviations of marks from curves are negligible. Note, we approximate the convection term in the admixture transport Eq. (14) with the use of the centraldifference scheme. This scheme gives a monotonic numerical solution of the convection-diffusion equation at the condition Peg  2 [41,46]; here, Peg is the local Peclet number based on mesh spacing. Otherwise, nonphysical oscillations can occur. In our terms, the maximal magnitude of Peg is estimated as Peg ¼ jwy jDy1 = 1.7 therefore the condition on Peg is satisfied. 5. Results of numerical simulations 5.1. A Homogeneous domain We briefly report on simulation of convective flows and mass transfer in a homogeneous isotropic porous domain considered

Fig. 7. The maximal fluid velocity jWjmax (a), relative admixture mass Rm (b), and depth of penetration hp (c) as a dependence on time in a domain without an interior layer and with an interior layer at properties (35)-(37).

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

as a basic case for comparing with convection in a domain containing an interior layer. The parameters are as declared in Eqs. (22), (25). One can roughly estimate the time of convection onset tth based on the threshold Rayleigh-Darcy number Rdth  40 in the concentration Rayleigh-Benard problem discussed in Section 4.2. In our problem, when an admixture starts to diffuse into a domain, the concentration boundary layer forms near the upper boundary. The width of this layer grows with time until reaches the threshold value hth estimated as hth ¼ Rdth ls D1 =ðagk1 qc Þ or, taking into  account Eq. (18), as hth ¼ Rdth H=Rd . We assume convection in a semi-infinite space analyzed here to start at the same magnitude of Rdth as in the case of bounded layer in the Rayleigh-Benard consideration. The time t th is related with hth according to the relation

791

The horizontal scale of downward flows increases whereas the boundary layer near the top starts after some period to generate new small fingers. To quantify an admixture propagation, some characteristics are exhibited in Fig. 7. As discussed, the module of maximal fluid velocity jWjmax starts to grow from about 150 days. Convection reaches a developed stage by about 400 days when jWjmax stops to grow sharply and starts to oscillate near a constant being equal about 40 mm d1. Convective flows near the upper boundary lead to increasing the admixture concentration gradients and,

2

t th ¼ hth =D1 that can be reduced to the following form used for esti2

2

mation: t th ¼ Rdth H2 =ðRd D1 Þ  550 days. We show below the module of maximal fluid velocity jWjmax (Fig. 7(a)) equal to zero if convection is absent and growing with convection onset. One can obtain from our simulation t th  150 days that is numerical and analytical threshold times are of the same order of magnitude. The numerical tth is shorter than analytical t th because our assumption of equal magnitudes of Rdth in the considered and RayleighBenard problems is not completely correct. Convection in a semiinfinite space has to start early at lower Rdth as there is no a bottom restricting convective flows. An initial stage of convection is shown in Fig. 5; a fragment of domain in an enlarged scale is presented. When the diffusion boundary layer near the top becomes unstable, the small brine fingers appear in random locations. They combine to a group with heading fingers in the center (Fig. 5(a) at x = 11  14 m) or form a series of nearly similar fingers located periodically (Fig. 5(b) at x = 14  20 m). Later, the behavior of fingers become strongly nonlinear. Convection becomes stochastic. The stage of developed convection is shown in Fig. 6. The irregular concentration pattern with the fingers moving downward and merging with each other in some places is observed. A convective behavior evolves in time.

Fig. 8. Field of admixture concentration in the domain with the interior layer at properties (35) at the times t = 280 (a), 730 (b), 1275 (c) days.

Fig. 9. Field of admixture concentration in the domain with the interior layer at properties (37) at the times t = 280 (a), 730 (b), 1275 (c) days.

Fig. 10. Field of admixture concentration in the domain with the interior layer at properties (36) at the times t = 280 (a), 600 (b), 730 (c), 1275 (d) days.

792

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

consequently, to enhancing mass transfer. The curve of relative admixture mass Rm shows that, in the beginning with Rm = 1, the admixture enters the domain only due to diffusion (Fig. 7(b)). Later, when convection becomes intensive, the amount of admixture in a domain enlarges several times comparing with that delivered by diffusion. The depth of admixture penetration hp exhibits a significant growth with time from about 350 days when convection approaches a developed stage (Fig. 7(c)). Density-driven convection in a homogeneous isotropic domain is discussed in detail in [6–8]. 5.2. Effect of properties of interior layer An interior layer located at the distance hd = 1.25 m to the top of domain is at lower porosity as the main medium. The width of layer hl = 0.25 m everywhere except Section 5.4 presenting a case of wider layer. We put three values of porosity /2 and, consequently, obtain three series of permeability k2 and effective coefficient of diffusion D2 calculated by Eqs. (8), (10). Appropriate dimensionless values are given below.

/2 ¼ 0:164;

j2 ¼ 0:505; r2 ¼ 0:485

ð35Þ

/2 ¼ 0:100;

j2 ¼ 0:099; r2 ¼ 0:398

ð36Þ

/2 ¼ 0:060;

j2 ¼ 0:020; r2 ¼ 0:325

ð37Þ

Eqs. (35)–(37) exhibit variations in /2 lead to strong changing in

j2 therefore the effect of different properties is associated mostly with the different permeability. The initial stage of fingering is the same as in Fig. 3 until convective flows reach the layer. Further development of convection at parameters (35) corresponding to light deviations from parameters of main medium is shown in Fig. 8. The interior layer is displayed in Figs. 8–14, 16 with gray color. The permeability of layer is lowered about twice relative to the main domain, however, it does not cause appreciable changing in the convection pattern. Movements of brine fingers downward are predominantly vertical. Some discrepancies with convection in a homogeneous domain are observed above the interior layer. Here, fingers slightly expand as their propagation throughout the layer is obstructed because of less permeability. The dynamical and mass transfer characteristics plotted in Fig. 7 are close to those in a homogeneous domain. Consider the case of layer at very different properties (37). The permeability is lowered in 50 times. Brine propagation throughout

Fig. 11. Fragments of fluid velocity (left) and concentration (right) fields in the case of interior layer at properties (35) at the times t = 280 (a), 730 (b), 1275 (c) days.

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

the layer is very difficult therefore, as shown in Fig. 9, an admixture is accumulated between the upper boundary and the layer filling this area densely. Only a part of admixture is able to travel down feeding thin fingers under the layer which grow in time. In the beginning of motion till the time about 300–400 days, the effect of interior layer is neglected: the maximal fluid velocity jWjmax in Fig. 7(a), admixture mass Rm in Fig. 7(b) and the depth of penetration hp in Fig. 7(c) are near those obtained in a homogeneous domain. After 400 days, the effect of interior layer becomes visible: the values in Fig. 7 decrease in a comparison with ones in the case of no layer. The interior layer appreciably constrains a propagation of convective flows reducing jWjmax and hp . Since an admixture is accumulated in the neighborhood of the upper boundary, the concentration gradients become smaller here. It leads to enhancing the mass transfer due to convection to a lower extent and decreasing Rm relative to the case of homogeneous domain. The case of moderate changing in layer properties corresponding to Eq. (36) is complicated and therefore considered in more detail. Concentration patterns are shown in Fig. 10. Brine flows transport an admixture down leaving its appreciable fraction above the interior layer. The concentration field between the layer and the upper boundary in Fig. 10(a)-(b) resembles the concentration field in Fig. 9(a)-(b) in a sense of dense filling this area with an

793

admixture. From Fig. 10(b), one can identify an ability of interior layer to pass flows throughout as weak till about 600 days. Later, the dynamic pattern changes crucially. Small fingers of pure water displayed with white color start to rise from the interior layer to the upper boundary. As exhibited in Fig. 10(c)-(d), the fingers of pure water split concentrated brine into a series of formations, which can be identified as roughly periodical ‘‘brine balls”. The average horizontal scale of such formations is about 3 m. Comparing patterns in Fig. 10(c)-(d), some distances between ‘‘brine balls” increase. As a result, some pure spots between the interior layer and the top not observed before (Fig. 10(b)) are visible later (Fig. 10(d)). Flows below the interior layer do not change dramatically and, as in the basic case (Fig. 6), are the falling predominantly vertical brine fingers which transport an admixture downward. The maximal velocity jWjmax is close to that in the case of very low permeability up to about 650 days (see curves at Eqs. (36) and (37) in Fig. 7(a)). Such behavior indicates a similar effect of interior layer on passing flows despite the values of permeability differ about five times. However, jWjmax at Eq. (36) increases sharply in a period from about 650 till 800 days and later exceeds the magnitude of jWjmax in the basic case of no layer. The extreme growth in jWjmax is associated with a decrease in porosity. As W ¼ U=/, where U is the filtration velocity to be continuous, the fluid

Fig. 12. Fragments of fluid velocity (left) and concentration (right) fields in the case of interior layer at properties (37) at the times t = 280 (a), 730 (b), 1275 (c) days.

794

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

velocity W increases by step if brine crosses the interface between media at high and low porosities. Analyzing the behavior of jWjmax , one can assume the regime of passing flows through the layer to change dramatically during the period between 650 and 800 days: the interior layer constraining effectively the motion of brine in early times does not do it later. The trend in behavior of Rm in Fig. 7(b) is in agreement with our assumption. After 800 days, the Rm curve deviates from one corresponding to the lowpermeable layer and approaches more and more the curve in the case of no layer. The depth of penetration hp shows the same trend as Rm (Fig. 7(c)). Dynamic and concentration patterns displaying in more detail a central (in x-direction) upper (in y-direction) fragment of domain are presented in Figs. 11–13. The fluid velocity vectors in the left figures are plotted in ten grid point intervals. As shown in Fig. 11 corresponding to the case of light changing in layer properties, convective motions in the layer are controlled by flows nearby during all time. A brine falls whereas pure water rises filling the space above the layer leaved by brine. This behavior is verified by a comparison of left and right figures: flows at downward velocity are associated with brine (dark color) whereas flows at upward velocity are attributed to water (light color). The dynamic behavior of brine in the interior layer is natural convection.

The low-permeable layer passes brine differently. As exhibited in Fig. 12, brine always moves here downward and no fluid is able to rise. As a result, there is no supply of pure water to the area above the layer. Note, the brine velocity in the layer is about vertical that follows from the boundary conditions at the interface (20) at dramatic changing in properties. If a brine flow travels from the main medium into the interior layer, the fluid velocity changes at the interface by step. The horizontal and vertical components are evaluated at parameters (37) as wx2 ¼ 0:067wx1 and wy2 ¼ 3:3wy1 . We see the horizontal component appreciably decreases while the vertical component increases. The decrease in wx2 is associated with the dramatically lowered permeability, however, the increase in wy2 is caused by the decrease in the porosity. On the contrary, the flow leaving the layer in the main medium spreads out along the interface just below the layer with the increased horizontal velocity. Released fluid under the layer feeds new brine fingers. Dynamic behavior of fluid in the interior layer is forced convection driven mostly by the pressure gradients rather than the gravity force. The dynamic and mass transfer patterns for the case of interior layer at moderate changes in properties are given in Fig. 13. Initially, the layer passes brine in a similar manner as a lowpermeable layer. As shown in Fig. 13(a) (left), the fluid motions

Fig. 13. Fragments of fluid velocity (left) and concentration (right) fields in the case of interior layer at properties (36) at the times t = 280 (a), 730 (b), 1275 (c) days.

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

795

the condition @ P=@y < 0 is satisfied at all coordinates x in all considered times. On the contrary, @ P=@y > 0 is observed if the properties of interior layer vary dramatically according with Eq. (37). The behavior of @ P=@y in the intermediate case of properties (36) is changeable. The values of @ P=@y are positive in the early times, however, later they turn to negative and become negative almost everywhere except some places. In those places, intensive downward flows traveling from ‘‘brine balls” occur. The typical behavior of pressure is illustrated by P profiles in Fig. 14; some fragments of P profiles along the central vertical are presented. As clear, inside the interior layer, the conditions @ P=@y < 0 at properties (35) and @ P=@y > 0 at (37) are satisfied in both times whereas the value of @ P=@y changes at (36) from positive at 500 days to negative at 1275 days. Consider the vertical component of Darcy’s Eq. (13). We have

  @P uy ¼ Rd j þS @y

Fig. 14. The vertical pressure profile in the center of domain with an internal layer at properties (35)-(37) at the times t = 500 (a), 1275 (b) days.

ð38Þ

For the density, the inequality S P 0 is always true. If @ P=@y > 0 everywhere, as it observed in the case of very different layer properties (37), the vertical velocity uy has to be negative: uy < 0. It means, that fluid in such layer must move downward everywhere. This is the case of forced convection. For natural gravity-driven convection to develop, the condition uy > 0 has to be satisfied in some points of layer. To obtain uy > 0, the inequality @ P=@y < S should be true. As S is positive, the last inequality can be reduced to the form @ P=@y < 0 being necessary but not sufficient. We find, the condition @ P=@y < 0 satisfied somewhere in the interior layer is necessary for natural convection to develop that is in agreement with data on @ P=@y in the cases of layer properties (35) and (36) discussed above. The transition from forced to natural mode of convection observed at (36) is associated with rearranging the pressure field in the interior layer. The curves of P in Fig. 14 are broken at the interfaces that is resulted from the Darcy’s law coupled with boundary conditions (20). 5.3. Effect of location of interior layer

have no upward direction in the area of layer. Later, some upward flows through the layer are arisen (Fig. 13(b) (left)). The time instant in Fig. 13(b) is in a period associated with transition from one convection regime to another predicted from the analysis of Fig. 7(a). We see that total downward fluid motions are rearranged into the complex motions directed both downward and upward depending on coordinates. If fluid moves in one direction, its velocity is restricted because of the conservation laws. When opposite direction arises, fluid is able to move with increased velocity as a space leaved by downward fluid is filled with rising fluid. So, transition from forced convection to natural convection characterized by the velocity increase is observed in a period from about 650 till 800 days. Fingers of pure water rising throughout the layer split brine in the area above into ‘‘brine balls”, as discussed in comments to Fig. 10. We see now, that formation of ‘‘brine balls” is associated with transition from the regime of forced convection to the regime of natural convection. One of ‘‘brine balls” is clearly shown in Fig. 13(c). A ‘‘brine ball” transformed into a finger below the layer provides an effective transport of admixture away from the upper boundary. A periodic structure of thermal formations resembling ‘‘brine balls” in shape is observed in steady Rayleigh-Darcy convection (thermal, not haline) in a porous domain divided by an extremely thin low-permeable layer [47]. To find conditions responsible for the mode of convection to take place in an interior layer, the pressure fields were examined. We focus on the vertical component of the pressure gradient @ P=@y because are interested in fluxes moving across the layer. We found that, in the case of light variations in properties (35)

The interior layer at properties (36) is considered. In the reference case discussed in detail early, the distance between the layer and upper boundary hd = 1.25 m. In this section, hd varies: hd = 0.75  5.0 m. As shown in Fig. 15(a), if the layer situates at hd = 0.75 m, the maximal velocity jWjmax decreases earlier than in the reference case because brine fingers reach the layer in a shorter time; the magnitude of jWjmax falls down to 25 mm d1 at about 250 days. As the velocity decreases, motions in the layer is forced convection. However, at about 450 days, the magnitude of jWjmax increases sharply nearly twice that indicates the transition from the forced to natural mode of convection in the layer. Mass transfer from the upper boundary characterized by Rm is lowered relative to the reference value, however, the discrepancy reduces with natural convection to develop (Fig. 15(b)). The behavior of depth of penetration hp is surprising: the curve of hp is close to one at hd = 1.25 m during an initial period and takes off after 650 days (Fig. 15(c)). The structure of convection looks like that in the reference case: between the top and interior layer, brine is split into ‘‘balls” transformed into long fingers below the layer (Fig. 16(a)). However, if the interior layer approaches the top, the characteristic scale of ‘‘brine balls” becomes shorter. The average horizontal scale of ‘‘balls” in Fig. 16(a) is about 1.6 m. If the interior layer distances from the top, its effect on convection is weakened. For example, at hd = 5.00 m, the behavior of jWjmax , Rm , and hp in Fig. 15 is close to that observed in a domain with no layer in Fig. 7. Perhaps, the distant layer will influent flows and mass transfer after a long time, however that far future is

796

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

Fig. 16. Field of admixture concentration in the domain with the interior layer at different values of distance to the upper boundary hd and width hl : hd = 0.75 m, hl = 0.25 m (a), hd = 5.00 m, hl = 0.25 m (b), hd = 1.25 m, hl = 0.80 m (c) at the time t = 1275 days.

Fig. 15. The maximal fluid velocity jWjmax (a), mass ratio Rm (b), and depth of penetration hp (c) as a dependence on time in a domain with an interior layer at different values of distance to the upper boundary hd and width hl .

beyond our consideration. As shown in Fig. 16(b), an admixture accumulates above the layer but the structure of ‘‘brine balls” does not appear.

5.4. Effect of width of interior layer A porous domain contains the wider interior layer comparing with that considered early. We put the layer width hl = 0.80 m and compare with the case of hl = 0.25 m. The properties of layer are defined by (36), the distance to the top hd = 1.25 m. As shown in Fig. 15(a), the maximal velocity jWjmax behaves like that in the reference case up to 650 days when the mode of convection into the layer starts to change. The transition between forced and natural convection in the wider layer lasts longer therefore the value of jWjmax increases slower than in the thinner layer. In the last case, after 800 days when natural convection develops completely, the

mass transfer is intensified while it does not in the first case. As observed in Fig. 15(b), the curve of Rm at hl = 0.25 m rises sharper as at hl = 0.80 m after 800 days. The depths of admixture penetration hp is lowered if flows travel throughout the wider layer, as it exhibited in Fig. 15(c). Fig. 16(c) demonstrates the convective structure formed by the time t = 1275 days which can be compared with Fig. 10(d) showing the appropriate pattern in the reference case. We see the wader layer allows the roughly periodic structure to form above itself as well. It is possible because arising upward flows start to transport pure water to the top. If the width of layer is bigger, the average horizontal scale of ‘‘balls” is longer. This scale is about 5 m in Fig. 16(c). Note, some ‘‘balls” feed not one but several fingers below the layer that is resulted from their longer scale.

6. Summary and conclusion Numerical simulations of stochastic haline convection in a porous medium in the context of admixture transport in geothermal systems were performed. Convection is driven by a denser brine fed by dissolved admixture, which diffuses from the upper boundary. The domain contains a horizontal interior layer at the lower porosity and related lower permeability; the coefficient of diffusion decreases into the layer as well. The objective of this work

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

was to study the effect of interior layer on convection dynamics and mass transfer in all domain. The effect is associated mostly with decreasing in the layer permeability. Two modes of convection depending on the magnitude of permeability are observed in the layer. If permeability is lowed lightly, fluid here is involved into natural gravity-driven convection with both falling and rising flows. In this case, the effect of layer is little: the admixture mass transported from the top is reduced very little whereas the maximal fluid velocity and depth of admixture penetration are about the same comparing with convection in a domain with no layer. Dynamic behavior of fluid in a layer at the dramatically lowed permeability is forced convection driven by the pressure gradients and consisting of only downward flows. This case is associated with a decrease in the fluid velocity, mass transfer and depth of penetration. In the case of moderate variation in the permeability, we have obtained an important result showing that the mode of convection in the layer can change in time. In an early period, fluid propagates throughout the layer by forced convection, which is transformed later into natural convection. The transition between first and second types of motion is characterized by increasing in the velocity in the layer and appearing upward flows, which transport pure water from the area below the layer to the area above it. Water supplied splits brine above the layer into roughly periodic ‘‘brine balls”. The mode of convection to occur is controlled by the vertical pressure gradients, which have to be positive everywhere for a development of forced convection and negative somewhere for a development of natural convection. Variations in a location of interior layer relative to the upper boundary show if the layer is close to the top, ‘‘brine balls” form earlier and have smaller scale. If the layer distances from the top, the pattern of ‘‘brine balls” does not observed, at least, in a considered period. An increase in the width of layer leads to slower transition between forced to natural convection and to a formation of longer ‘‘brine balls”. Our study demonstrates the interior layer is able to order convective flows above itself and is responsible for the scale of roughly periodic structure formed. Analyzing this structure, it is possible to obtain some information on the layer location and properties. Conflicts of interest The author declares that there are no known conflicts of interest. Acknowledgements The author would like to thank G.G. Tsypkin for many fruitful discussions. This work has been supported by the Russian Science Foundation, Russian Federation (grant № 16-11-10195). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.ijheatmasstransfer. 2018.08.019. References [1] D.A. Nield, A. Bejan, Convection in Porous Media, Springer, 2006. [2] J. Bear, A. Cheng, Modeling Groundwater Flow and Contaminant Transport, Springer, 2010, 834 p. [3] L. Anderson, New Theory of the Earth, Cambridge University Press, Cambridge, 2010. [4] E. Soboleva, Numerical simulation of haline convection in geothermal reservoirs, J. Phys. Conf. Ser. 891 (2017) 012105.

797

[5] B. Metz, O. Davidson, H. Coninck, M. Loos, L. Meyer, IPCC Special Report on Carbon Dioxide Capture and Storage, Cambridge University Press, Cambridge, 2005. [6] A. Riaz, M. Hesse, H.A. Tchelepi, F.M. Orr, Onset of convection in a gravitationally unstable diffusive boundary layer in porous media, J. Fluid Mech. 548 (2006) 87–111. [7] R. Farajzadeh, H. Samili, P.L.J. Zitha, H. Bruining, Numerical simulation of density-driven natural convection in porous media with application for CO2 injection projects, Int. J. Heat Mass Transf. 50 (2007) 5054–5064. [8] M. Bestehorn, A. Firoozabadi, Effect of fluctuations on the onset of densitydriven convection in porous media, Phys. Fluids 24 (2012) 114102. [9] H.E. Huppert, J.A. Neufeld, The fluid mechanisms of carbon sequestration, Annu. Rev. Fluid Mech. 46 (2014) 255–272. [10] T.F. Faisal, S. Chevalier, Y. Bernabe, R. Juanes, M. Sassi, Quantitative and qualitative study of density driven CO2 mass transfer in a vertical Hele-Shaw cell, Int. J. Heat Mass Transf. 81 (2015) 901–914. [11] M. Paoli, F. Zonta, A. Soldati, Dissolution in anisotropic porous media: Modeling convection regimes from onset to shutdown, Phys. Fluids 29 (2017) 026601. [12] A.W. Islam, M.A.R. Sharif, E.S. Carlson, Numerical investigation of double diffusive natural convection of CO2 in a brine saturated geothermal reservoir, Geothermics 48 (2013) 101–111. [13] A. Islam, A.K.N. Korrani, K. Sepehrnoori, T. Patzek, Effects of geochemical reaction on double diffusive natural convection of CO2 in brine saturated geothermal reservoir, Int. J. Heat Mass Transf. 77 (2014) 519–528. [14] C.M. Oldenburg, K. Pruess, Layered thermohaline convection in hypersaline geothermal systems, Transp. Porous Media 33 (1998) 29–63. [15] R. McKibbin, Modeling heat and mass transport processes in geothermal systems, in: K. Vafai (Ed.), Handbook of Porous Media, Taylor & Francis Group, Boca Raton, 2005, pp. 545–571. [16] H. Emami-Meybodi, Dispersion-driven instability of mixed convective flow in porous media, Phys. Fluids 29 (2017) 094102. [17] M.Z. Saghir, M.R. Islam, Double diffusive convection in dual-permeability, dual-porosity porous media, Int. J. Heat Mass Transf. 42 (3) (1999) 437–454. [18] A.H.-D. Cheng, Multilayered Aquifer Systems: Fundamentals and Applications, Marcel Dekker Inc., New York, 2000, p. 400. [19] S. Rapaka, R.J. Rawar, Ph.H. Stauffer, D. Zhang, Sh. Chen, Onset of convection over a transient base-state in anisotropic and layered porous media, J. Fluid Mech. 641 (2009) 227–244. [20] D. Daniel, A. Riaz, H.A. Tchelepi, Onset of natural convection in layered aquifer, J. Fluid Mech. 767 (2015) 763–781. [21] V.A. Shargatov, Instability of a liquid-vapor phase transition front in inhomogeneous wettable porous media, Fluid Dyn. 52 (1) (2017) 146–157. [22] G.G. Tsypkin, Stability of the evaporation and condensation surfaces in a porous medium, Fluid Dyn. 52 (6) (2017) 777–785. [23] P. Xu, B. Yu, Developing a new form of permeability and Kozeny-Carman constant for homogeneous porous media by means of fractal geometry, Adv. Water Resour. 31 (2008) 74–81. [24] M. Matyka, A. Khalili, Z. Koza, Tortuosity-porosity relation in porous media flow, Phys. Rev. E 78 (2008) 026306. [25] A.A. Martins, P.E. Laranjeira, C.H. Braga, T.M. Mata, Modeling of Transport Phenomena in Porous Media Using Network Models, in: Kong Shuo Tian, HeJing Shu (Eds.), Progress in Porous Media Research, Nova Science Publisher, Columbus, Ohio, 2009, pp. 165–253. [26] P. Guo, Dependency of tortuosity and permeability of porous media on directional distribution of pore void, Transp. Porous Media 95 (2012) 285–303. [27] J.M.P.Q. Delgado, A critical review of dispersion in packed beds, Heat Mass Transf. 42 (2006) 279–310. [28] J.L. Musuuza, F.A. Radu, S. Attinger, The effect of dispersion on the stability of density-driven flows in saturated homogeneous porous media, Adv. Water Resour. 34 (2011) 417–432. [29] J.L. Musuuza, F.A. Radu, S. Attinger, The stability of density-driven flows in saturated heterogeneous porous media, Adv. Water Resour. 34 (2011) 1464– 1482. [30] Z. Jamshidzadeh, F.T.-C. Tsai, S.A. Mirbaghery, H. Ghasemzadeh, Fluid dispersion effects on density-driven thermohaline flow and transport in porous media, Adv. Water Resour. 61 (2013) 12–28. [31] L. Wang, Yu. Nakanishi, A. Hyodo, T. Suekane, Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure, Int. J. Greenhouse Gas Control 53 (2016) 274–283. [32] G.B. Rossa, K.A. Cliffe, H. Power, Effects of hydrodynamic dispersion on the stability of buoyancy-driven porous media convection in the presence of first order chemical reaction, J. Eng. Math. 103 (2017) 55–76. [33] J.-H. Kim, A. Ochoa, S. Whitaker, Diffusion in anisotropic media, Transp. Porous Media 2 (1987) 327–356. [34] H. Chou, L. Wu, L. Zeng, A. Chang, Evaluation of solute diffusion tortuosity factor models for variously saturated soils, Water Resour. Res. 48 (2012) 10539. [35] R.J. Millington, J.P. Quirk, Permeability of porous solids, Trans. Faraday Soc. 57 (1961) 1200–1206. [36] E.B. Soboleva, G.G. Tsypkin, Numerical simulation of convective flows in a soil during evaporation of water containing a dissolved admixture, Fluid Dyn. 49 (5) (2014) 634–644. [37] E.B. Soboleva, G.G. Tsypkin, Regimes of haline convection during the evaporation of groundwater containing a dissolved admixture, Fluid Dyn. 51 (3) (2016) 364–371.

798

E.B. Soboleva / International Journal of Heat and Mass Transfer 127 (2018) 784–798

[38] S. Patankar, Numerical Heat Transfer and Fluid Flow, Publishing Corporation, New York, 1980. [39] Kh. Aziz, A. Settari, Petroleum Reservoir Simulation, Applied Sciences Publishers Ltd., London, 1979. [40] A.A. Alexandrov, K.A. Orlov, V.F. Ochkov, Thermophysical Properties of Thermal Power Engineering Working Substances: reference book, MPEI Press, Moscow, 2009. [41] D.R. Lide (Ed.), CRC Handbook of Chemistry and Physics, Taylor and Francis, Boca Raton, 2006. [42] A.E. Gill, Atmosphere-Ocean Dynamics, Academic Press, San Diego, 1982. [43] C.A.J. Fletcher, Computational Galerkin Methods, Springer-Verlag, New York, Springer series in computational physics, 1984.

[44] L.D. Landau, E.M. Lifshitz, Fluid Mechanics, second ed., Pergamon, New York, 1987. [45] G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review, McGraw-Hill 1130 (1968). [46] J.H. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Springer, Berlin, 2002, p. 423. [47] D.R. Hewitt, J.A. Neufeld, J.R. Lister, High Rayleigh number convection in a porous medium containing a thin low-permeability layer, J. Fluid Mech. 756 (2014) 844–869.