Seasonal heat storage in the soil

(1)

(2)

n

ISBN 90 6231 147 4

DUTCH EFFICIENCY BUREAU - PIJNACKER

1985

ASQMAL UE/AT STORAGE DS

im

)L

PROEFSCHRIFT

ter verkrijging van de graad van doctor in de

technische wetenschappen aan de Technische

Hogeschool Delft, op gezag van de rector

magnificus prof.dr. J.M. Dirken, in het

openbaar te verdedigen ten overstaan

van het College van Dekanen op

donderdag 19 december te 14.00 uur

door

Gerardus Antonius Maria van Meurs

wiskundig ingenieur

geboren te Wateringen

TR diss

1468

(3)

Di •; proe-f schri-f t is goedgekeurd door de promotor Pro-f . ir. C.J. Hoogendoorn

aan mijn ouders voor Ada

(4)

-CHAPTER 1

©EMElRiAL IMimODUCTDOM

1.1 Introduct:in

Man was able to heat up his shelter with wood when he discovered the principle of making fire. At the same time he must have become aware of the fact that collecting this wood was less strenuous and dangerous in summertime (or daytime) than in wintertime (or nightime) because of weather conditions and road conditions, so he initiated the method of energy storage. Wood was replaced by fossil fuels (coal, oil and natural gas) perhaps a century ago and since people changed over from wood to these fossil fuels it has been possible to store the energy just by storing the the material in which the energy is present. In fact, the energy has been stored in a concentrated form undernearth the earth mantle for centuries in this way.

The interest in energy storage directly as sensible or latent heat, however, is rather new. It really started after the first oil crisis of 1973-1974. People in the western countries then not only became aware of their dependence on the resources of fossil fuels, but they also realised that these natural resources are finite. This concept was already stressed by the Club of Rome ("Limits to Growth"). It became apparent that these stored, concentrated energy forms are being used at such a rate that they will be depleted in the not-too-distant future. Furthermore, during the last decade, the cost of conventional fuels has increased rapidly and also there have been political uncertainties related to the supply of oil. As a result a growing interest occured in new energy resources. The most significant of these are nuclear and solar energy. The main drawback of nuclear energy is that not all the problems are solved in managing its radioactive waste. Solar energy, on the other hand, is available without causing any pollution. It is transmitted from the sun by electromagnetic radiation. Solar radiation can be converted into electricity by using solar

(7)

12

-cells and into heat by application of solar collectors. A disadvantage o-f solar radiation is that it varies during the day and with the season. Also the demand -for- energy varies. However, the peak values for both do not coincide with each other. Useful solar collection is possible at times when there is little or no demand for energy. On the other hand, there are times when useful solar energy can not meet the energy demand. A phase shift

occurs and energy storage might come in useful .

Heat storage can be applied whenever the moment of production and demand of heat are out of phase. These phase shifts are significant both for daily and seasonal time scales for solar heating systems.

Short term (ST) storage is usually sized such that the system is able to cover the heat load within a period of one or two days. For instance, heat is collected at day time and used for the night. It is a conventional technique in many applications. Typical examples of ST storage of sensible heat (storage in a material by virtue of its temperature rise) are tanks and basins filled with water, e.g. a water tank connected to a solar boiler for domestic hot water supply. Storage can also be realised by using the latent heat of materials undergoing a phase change (e.g. paraffin wax). Such materials can store more energy per unit volume than sensible heat storage materials. However, they have some predominant disadvantages (undercooling, segregation, low thermal conductivity in the solid phase).

Long term (LT) storage or seasonal storage covers the effect of seasonal meteorological variations on the performance of a solar space heating system. When such a storage is added to a solar energy system, the next question has to be answered during the design of the system. How large must the collecting array and the storage volume be? Technically, it will always be possible to have such a large collecting array and storage volume to carry the system over all periods of inadequate supply of solar energy. The system will then be so large that it is not economical 1y viable at all. The alternative is to accomplish only a percentage of the energy demand and to have an auxiliary energy supply to replenish the deficit.

Adding this storage concept to a solar heating system is not the only application of LT heat storage. Besides solar heat, examples of heat production which might be used for storage are heat supplied by power stations or waste heat supplied by industry. In order to reduce the heat loss per stored energy unit, it is preferable to make one big store instead

13

-of a number -of smaller stores together having the same volume. The consumption of the stored heat can take place by the heating demand of dwellings or greenhouses. If the heat is destined for dwellings then the storage concept must be integrated in a district heating system.

1.2 Heat storage in the soil

This thesis reflects an investigation on the thermal performance of the soil as seasonal heat storage medium. The thermal characteristics of the soil storage and the heat transfer from a heat exchanger in the soil will be discussed. This investigation forms a contribution to a project to realise such a storage in a solar heating system in the Netherlands.

Meanwhile such a solar system is realized. This system serves as a district heating system for a new extension of the town of Groningen, called Beyum, which is situated in the north-eastern part of the Netherlands. The project includes 96 dwellings. The solar energy system delivers a substantial amount of the yearly demand for space heating and for domestic hot water. It is expected that approximately 65 'A of the yearly heating demand of the

dwellings will be met by the solar energy system. For space heating there is a central auxiliary heater, whereas for the domestic hot water provision each dwelling has its own back-up system. Each dwelling is fitted up with 25 m2 tubular collectors, in which a vacuum prevails, on its roof. The dwellings are connected to a central storage reservoir by a distribution network. In summer time the surplus of heat coming from the solar col 1ec tors is flowing through this network towards the reservoir. In wintertime the heat withdrawn from the reservoir to meet the heat demand of the dwellings is flowing towards the dwellings. The reservoir is a cylindrical part of the soil which acts as the storage volume for sensible heat. The volume of the reservoir is approximately 23980. m3 with a thermal capacity of 6.ó 161 0

■J/K. A heat exchanger is brought into the reservoir to transfer the heat to and from the soil. An insulating layer is placed above the reservoir to suppress heat loss to the atmosphere. At the other boundaries of the reservoir the soil itself has to be the insulator.

At most places in the Netherlands the phreatic face is just below ground level (8.5 till 2. m ) . As a result the soil in such a reservoir will be saturated with water. Density differences of the groundwater will occur when heat is introduced into the soil. Natural convection (free convection) then might develop. Also regional groundwater flow (forced convection)

(8)

-

14

-might take place when there are hor izon tal potent i al di fferences because of drainage, groundwater extraction and/or infiltration of rain. This is the reason that diffusion as well as convection of heat are responsible for a transfer of heat from the reservoir to the surroundings. Both these phenomena will then influence the heat loss.

t.3 Outline of this thesis

The various concepts of seasonal heat storage will be outlined in Chapter 2. After that, the theory of natural convection in porous media viill be di scussed in Chapter 3. Literature about this topic will also be presented in this chapter.

The governing equations describing the combined fluid floe and heat transpor t in porous medi a together with initial and boundary condi tions result in the mathematical model of our problem. The theory behind these equations and the mathematical model itself will be presented in Chapter 4. This mathematical model has to be solved in order to quantify several

aspects of the seasonal heat storage concept. It has been solved with a numerical technique. The applied numerical technique is a finite difference formu1 at ion and it wi11 be ment ioned in Chapter 5.

Computational results for a general case of seasonal heat storage will be presented in Chapter 6. For the description of the viscosity and the

density of water as a function of temperature two possibilities for each property have been considered. The temperature dependency of the viscosity is neglected and a constant value is taken or the temperature dependent viscosity is represented by tabulated values. The temperature dependency of the density can be linearized by defining a thermal expansion coefficient, which is constant, or the temperature dependent density can again be represented by tabulated values. The effect of these descriptions of the viscosity and the density of water on the development of natural convection will be discussed in this chapter. The influence of the permeability of the soil on this development and thus on the heat loss from the reservoir is investigated. A solution to suppress natural convection will also be discussed in this chapter. The influence of several system parameters on the thermal performance of the reservoir will be presented. Finally, the effect of a regional groundwater movemen t wi11 be an objec tive in thi s chapter.

The problem of the heat transfer between a vertic' tube of the heat

15

-exchanger and water-saturated soi 1 , cal 1ed the 1 ocal process of heat transfer, is illuminated in Chapter 7. A special developed numerical grid has been applied in order to describe the mechanism of heat transfer accurately. Stationary as we 11 as transient thermal behaviour of the transfer fluid are considered. The effect of convection on the heat exchange for the stationary case has been investigated. A ccmparison between the numerical results and an analytical approximation is given. After that, the total heat exchange within the reservoir is described. The conditions within these calculations for this investigation are chosen in agreement with the situation at the Groningen Project. At the end of this chapter the consequences for the description of the exchanging capability of the soil heat exchanger of the Groningen Project are presented.

First of all, a detailed description of the Groningen Project will be presented in Chapter 8. The description of the thermal behaviour in the reservoir and its surroundings, called the global process of heat transfer, together with a formulation for the heat exchange in the reservoir based on the results of Chapter 7 can now be given. Several new concepts to store the heat in the reservoir more efficiently are also considered. Each of these concepts can be realised by interconnecting the strings of the soil heat exchanger in a different manner. Their respective effects on the accomplished heat demand and the ratio of the injected and withdrawn heat are given.

A compar i son between measured data and compu tat iona 1 results on the Groningen Projee t will be given in Chapter 9. Also the results will be presented of an experiment carried out in the laboratory to investigate the 1ocal heat exchange between a vertical tube and the surrounding water-saturated sand.

(9)

17

-CHAPTER 2

SEASONAL HEAT STORAGE

2.1 Introduction

In most s o l a r energy a p p l i c a t i o n s di-fferences e x i s t between the time d i s t r i b u t i o n o-f a v a i l a b l e s o l a r energy and of energy demand. For space h e a t i n g systems, these d i - f i e r e n c e s are s i g n i f i c a n t on both d a i l y and seasonal time s c a l e s . Short term s t o r a g e i s u s u a l l y s i z e d such t h a t the system i s able to meet the load w i t h i n the p e r i o d of one or a few days. Long-term (LT> s t o r a g e or seasonal s t o r a g e , on the other hand covers the e f f e c t s of seasonal me t e or-o 1 ogi cal var i at ions on the performance of the s o l a r system.

Options of LT thermal s t o r a g e a r e :

- Water tanks and s o l a r ponds - Rocks

- 3oil storage - Aquifers

These options were discussed by Duffie 3; Beekman (I960) and Torrenti (1988)* Givoni (1977) gave an outline of the technical aspects of these options with their merits and objections. A study of the cost comparison of some of these options can be found in Nichols (1978). Svec & Palmer- (1980) and Torrenti (1930) remarked that solutions to use the undisturbed soil for LT-storage are economically the most promising ones. A general outline of subsurface heat storage can be found in Margen (1931) and Berntsson (1983). Eriksson (1933) gave a summary of the social and national aspects concerned with the introduction of storage systems for space heating. In this Chapter the first two options will be presented only shortly, the others will get more attention.

(10)

18

-2.2 Water tanks and solar ponds

Lunde (197?) investigated the per-formance of a solar heating system when an annual storage was utilized. The storage applied was a water tank. The thermal capacity was 3.2 M.J/K -for a square meter collector surface. To describe the all day heat collection, he integrated the collector equation o-f Hottel & Whillier, paragraph 4.3*3.1* for a period of time. He presumed that the storage tank is well mixed, that the collector inlet temperature equalled the storage temperature and that this temperature rose or fell linear with time. The resulting equations accurately predicted in a single step the per-formance of the annual storage system -for an entire month. The difference in predicted annual performance was only two percent, rel alive to an hourly simulation.

Tabor i 1980) described the storage capability of solar ponds. A solar pond

is a mass of water used as a solar collector. This mass of water is relatively shallow. The solar irradiation is absorbed at the bottom. In a normal pond, the heated water would rise to the top and the heat would be dissipated to the atmosphere. By imposing a density gradient, by dissolution of salt, the concentration increasing with depth, convect ion might be suppr-essed.

Braun et al (1981) developed a mathematical model to investigate the effect of the presence of a water tank as a seasonal storage of energy in solar-heating. The rate o-f useful energy gain of the collector is determined by the Hottel-Whil1ïer equation. The instantaneous energy balance of the building is considered as a first order process with a lumped effective building capacitance. The instantaneous heat loss rate of one house is considered to be proportional to the temperature difference between the insi de and outside. The propor tionali ty factor is assumed to be cons tan t. The same description is used for the rate of the heat of each compartment of the storage which is fully mixed. Their study has proved that stratificatvon effects are not important for seasonal storage space heating systems. For a water tank with forty compartments each with a uniform temperature the improvement of the annual solar fraction is less than 27.. They have also found that for a fixed annual fraction of the load met by solar energy significant reductions in the collector area can be realized by adding storage capac i ty. Greater reductions in col 1ec tor area requ iremen ts with increasing storage occur for higher solar fractions. Storage capacity

19

-and collector area are interchangeable in achieving an annual load met by solar energy. Then the cost-effectiveness will determine the ultimate choice.

2.3 Rocks

Large excavations in rock ■'caverns and pits) can be used for heat storage. In most cases these volumes are filled with water. Sometimes, these volumes can only be excavated partially. A block filled system is then created (Bjurstrom, 1933). Rocks can also be used as a storage medium when a great number of boreholes are filled with plastic tubes in which water flows <Andersson 't<_ Eriksson, 1939;. Also a combination of a cavern and boreholes

is possible ïNilsson, 1933). Then the cavern may act as a short term storage. The affection of the chemical alteration of the water composition in caverns on the efficiency can be found in Claesson & Ronge (1933). If this water is pumped through a heat exchanger then precipitation of earlier dissolved species might occur and the efficiency decreases.

2.4 Soi1 storage

When using the undisturbed soil as a storage medium of sensible heat, the heat exchange takes place by means ot a heat exchanger. This exchanger consists of vertical tubes or strings ''.U-tubes) which are inserted into the soil from the surface. Above the soil they are interconnected.

Shelton (1975) studied the heat loss by conduction of an underground storage of hemispherical geometry. He used an analytic steady state solution to determine the daily heat loss which was only 0.86154 V, of the storage

capacity. The thermal conductivity of the soil was 2. W/m.K. The radius of the storage region was 16. m and there was a temperature difference of 40*C between the storage and the undisturbed soil. The time required for the average heat-loss rate of a new system (starting with the surrounding soil with its natural temperature distribution) to approach the steady state value was approximately a year. He also concluded that the soil surrounding the heat storage region contributed y^ry little to the system's storage

capacity neither on a daily nor on a seasonal time scale.

Hustacchi et al (1931) developed a procedure for simulation and evaluation of seasonal storage of solar energy. A set of pipes inserted into their storage region took care of the heat exchange. They described the whole

(11)

28

-system by means o-f the behaviour of a cylindrical region around each pipe. They formulated the following hypotheses i) only conduction of heat takes place ii> the cylindrical region does not exchange heat with the surrounding soil. They found that the maximum possible heat fluxes from a vertical tube to the soil were obtained when the tubes are large and located as far away from each other as possible.

Platell (19811 described the addition of a ground heat storage to a solar heating system for which the maximum temperature of the store is not higher than 35. ""C (Sunstore Project). Then it is possible to use simple solar-collectors which are integrated in the design of the building. The heat storage is created by constructing vertical channels for the circulation of the water. The spacing depends on the type of ground and the depth depends on the size of the object together with technical restrictions. His 1 «4 temperature philosophy requires a concept of integrating the heating surface into the floor and/or the ceiling of the building in order to obtain sufficient area for the heat exhange with the room.

Cl aesson in Hellstrom (1931) made a model study of duct storage systems.

They used the term "duct" for a buried layer of pipes, a bundle of vertical pipes in clay, a system of boreholes in rock or a row of plate heat exchangers. They arranged a coupling between the local heat transfer process around a duct and the total heat transfer mechanism of the storage region and the surrounding ground in their numerical model. The local heat transfer mechanism includes the thermal characteristic of the heat transfer-fluid within the duct. Only heat conduction is taken into account for this study. Besides the energy efficiency or energy recovery factor ^p (ratio between the withdrawn heat and the heat input) they introduce the temperature efficiency as:

f - T

out 6

1| = (2.1) in 8

They computed the energy and temperature efficiency of a reservoir with radius R and height H after several cycles of injection and extraction of heat. The desired magnitude of" heat flux to and from the reservoir is 1 MW at a temperature level T, of the heat transfer fluid which has to be within the range 58 < T , < 98'C tx = 1 W/m.K and <Pc) ~ 2.18& J/mS.K].

21

-Only the results after the fifth cycle of varying the dimension of the

Table 2.1 Energy and temperature efficiency.

1

1/8 1 8 1 1 H 2 5 . 5 8 . 188. 31.5 79.4 R 12.5 2 5 . 5 8 . 31 .5 19.85 E B.47 0.75 8.83 0 .73 8.73 1 T 8.72 8.74 8.76 8.74 8.74

reservoir are given. The volume U is normalized on the case H=59. m and R=25. m. Their main conclusion is that the size of the storage is of decisive importance to the heat loss and storage efficiency. The heat loss compared to the storage capacity becomes excessive for small systems. Hellstrcm (1983) mentioned that besides the local and global heat transfer-mechanisms also a certain "steady flux" solution has to be taken into account. This contribution represents the slew redistribution of heat between various regions of the storage due to the circulation of the transfer fluid through these regions of different temperature levels. He also mentions the simulation of various small duct storage systems carried out in Sweden. The agreement with measurements is quite satisfactory. Olsson < 1933) mentioned a seasonal heat storage project in Sweden (Kullavik Project) for which the maximum temperature of the storage is óe.^C. The soil at the site of the storage consists of clay. The heat storage is divided into a region in the centre for storing heat at a high temperature surrounded by a region in which heat of a 1ower temperature will be stored. Heat withdrawn from the high temperature region (HTR) can be used directly for heating purposes while heat withdrawn from the low temperature region (LTR) will be raised to a higher temperature by means of a heat-pump. Consequently, radial heat loss from the inner region will be withdrawn in the LTR. The heat exchanger between the fluid and the clay consists of poly-ethylene pipes which are put into the ground vertically. The total length of tubes is 5696. m. The outer diameter of the tubes is 32. mm. Olsson (1934) added that the HTR contains 138 U-tubes which are placed in a square with a distance of 8.5 m between them. In the LTR the tubes are placed into three parallel lines with a distance of 2. m between them, while

(12)

22

-the distance between -the tubes is 1.5 m within -the lines. The distance between the outer tubes o-f the HTR and the inner line o-f the LTR is 5. m. The depth of the HTR is 3. m and the LTR is 12. m deep. The total storage volume is 8188. m5 of which the HTR takes up 280. ms. It operates as a short term storage during the summer. The number of dwellings connected to this system amounts to 53. The col 1ec tors are integrated in the construction of the roof and they have a total area of 540. m£ which are facing south east. They have a corrugated acrylic-plate as a cover. Two third of the net-heat demand will be met by solar energy while electrical energy for the heat pump will contribute by a fourth. The deficit will be supplied by oil.

In Holland Uan Koppen et al (1973) were the first who mentioned the addition of a seasonal storage in the ground to a solar heating system. They particularly went into the subject of thermal stratification in the storage reservoir. They expected that a systematic use of it generally leads to an increase in the heat gain of solar heating installations of approximately 15/1. In Fischer et al C197?) a first order engineering model was formulated to desribe the system. All unsteady heat transfer phenomena were reduced to quasi-stationary problems in this model. The capacity of the storage was formulated proportional to the difference between the upper and the 1ower limit of the temperature in the core of the storage. Their proportionality factor consisted of the heat capacity of the soil and the volume of the storage which was represented by the cubic of a characteristic length. They concluded that increasing the storage temperature, in general, leads to a decrease of the energy efficiency when the capacity of the storage is kept constant. However, when the temperature difference is small, this might be reversed.

Wijsman te Den Ouden (1936.) developed a model to simulate the heat flows in a

solar heating system for a group of solar houses. This model consists of three parts. The first one simulates the solar heating system of one dwelling. The second one simulates the heat flow and the heat loss in the distribution network. Finally, the last one simulates the performance of the seasonal heat storage. This part is a one-dimensional simplified model. The soil is represented by several shells between which only diffusive heat exchange takes pi ace. The heat f 1ow exchanged with the reservoir is

23 -represented by the next f uric t ion ,

q = « ♦ <P c ) (7 -T ) (2.2) p f r , in H

in which Tn is the soil temperature in the centre of the reservoir, 7 ■ is the inlet temperature of the fluid and i is the effectiveness of

the soil heat exchanger. This heat flow is distributed uniformly over the volume of the reservoir. When heat is withdrawn, this will only be done within the shells of the reservoir of which the temperature is above T^ . . The effect of a variation of the capability of several components of the solar system on the yearly yield of the system have been studied with their model. They concluded that for low temperature heating systems, the addition of a seasonal heat storage to a solar heating system results in a significant increase in performance only when so called high performance solar col 1ec tors are used, the distribution network ls very wel 1 insulated

and the distribution network has a length as smal 1 as possible. Increasing the number of dwellings connected to the reservoir above 286 in order to reduce the total heat loss in the system is not effected due to the fact that for this number of dwellings the heat loss of the reservoir is already a factor 3 smaller than the heat loss in the distribution network. Some of the results of this feasibility study can also be found in Wijsman < 1930) . This study has been done for the Groningen Project. It is found that the solar contribution from the LT storage increases strongly with increasing collector area.

Results of simulating the test period of the Groningen Project with this model can be found in Wijsman bc Uan Meurs (1984). After a modification of

the function which simulates the heat withdrawal, the model could follow the changes of the heat flow that is withdrawn. These changes were caused by a variation of the inlet temperature of the fluid to the soil heat exchanger. The modification is that the withdrawn heat flow depends an the temperature near the edge instead of the temperature in the centre of the reservoir.

In all the foregoing studies about heat storage in the soil diffusion is included and convection is 1ef t out of account. In the Nether 1ands the graundwater table comes at most places to about 1. m below the ground level. The soil in such a reservoir is then always fully saturated with water and together with conduction also convection in the groundwater is responsible

(13)

24

-for the transport of heat from the storage region to the surrounding soil.

2.5 Aquifers

An aquifer is best defined as a saturated permeable geologic unit that can transmit significant quantities of water under ordinary hydraulic gradients (Freeze it Cherry, 1979). This layer may be enclosed by layers of low

permeability (e.g. clay). Then it is called "confined".

The interest in the concept of heat storage i n an aqu i fer has increased during the last ten years. The basic idea is to inject hot water into an aquifer during periods of low energy demand, and then, when energy demand is high, to extract the water and to use the heat. This concept of heat storage will get more attention here than the first two options because natural convection as well as forced convection are major heat transport phenomena in this application of heat storage. The mathematical description of these phenomena is similar to the one used for our problem of heat storage in the upper soil which is water-saturated.

T ^ T+AT ^supply well unconfined aquifer confining layer/ ground level injecfion well

Figure 2.1: The concept of a doublet well system

Uarious examples of heat production can be used for heat storage in an aquifer. Besides solar energy, it can be waste heat coming from industry or a power station. An aquifer naturally provides a large volume of water as the storage medium at relatively low costs. Water is extracted from the aquifer through a supply well, and after being heated, it is injected into an injection well in the same' aquifer (e.g. figure 2.1). After the hot water has been produced from the injection well and the heat is used, the water will be reinjected into the supply well.

25

-6yuk (1983) presented an outline of aquifer thermal energy storage in the United States. He mentions two experiments of the injection of hot water (Mobile in Alabama and St Paul in Minnesota). The experiments on the first site are carried out by the Auburn University. Initially, when using a single well system, the test in Alabama was not very successful since massive clogging of the formation due to swelling and dispersion of clays developed upon injection. Incompatibility of water from two aquifers seemed to be the cause. The concept of a well doublet was then realised to solve this problem. The horizontal permeability of the aquifer is 6.2 I0"1 1 m£. The horizontal to vertical permeability ratio (anisotropy ratio) amounts to 6.7. For the first cycle, the average temperature of the injected water was 5S.5*C. The amount of injection was 2Ó009. m3 and the rate of injection was 45. m3/hr (Melville et al , 1981). After one month of storage, recovery of equal volume yielded an energy recovery factor of 0.55. Tracer studies indicated small di spersi on but significant natural convect i on due to buoyancy forces caused by the density differences of the water. The temperature drop during recovery was almost linear. During the second cycle, the average temperature of injection was raised to 81.*C. The amount of injection was 53000. mE of hot water. The temperature drop during recovery was again linear. However, it became clear that continued withdrawal would lead to a final energy recovery factor of only 0.46. Apparently, buoyancy had increased and more relatively cold water was reaching the well. When only withdrawal from the upper half of the well

took pi ace, the temperature jumped by 3. *C. Subsequent development indicated a recovery factor of perhaps 6.47.

At the project in St Paul (Minnesota), superheated water is injected into deep aquifers some 200. meters down. Since temperatures are around 156."'C,

the water will be under pressure.

The concept of aquifer storage of thermal energy has been studied extensively at the Lawrence Berkeley Laboratory since 197Ó (Tsang et al, 1974, 1973, 1981). For the simulation of the thermal performance around the injection well of a doublet system, radial symmetry with respect to the injection well was considered. This is justified if the spacing between the wells is large enough so that the thermal behaviour around this injection well is not significantly affected by the neighbouring well. Detailed simulation of the experiment in Mobile (Alabama) has been carried out to

(14)

demonstrate the predictive ability of the numerical model. With this model several schemes are bei nq expi ored to opt imi ze the energy recovery. These schemes include a partial penetration in the lower half of the well during injection to reduce tilting of the thermal front and a partial penetration in the upper half of the well during production to extract only hot water which resides in the upper part of the aquifer due to natural convection. The effect of natural convection can be neglected on the one hand in aquifers with a 1ow permeability or on the other hand when a short length of the storage cycle or when a smal1 temperature difference between the water injected and the ambient water is present. Several parameter sensitivity calculations were made. It was found that although the recovery factor can be made to reproduce the experimental data, for many alternatives the time rate of the decrease of the calculated production temperature is always faster than the experimental value. Their preliminary conclusion is that this discrepancy may be due to either thermal dispersion or aquifer heterogeneity. Due to the effect of changing the strategy for the production period of the second cycle of the experiment at Mobile (Alabama), it seemed useful to study alternative injection and production schemes to maximize the recovery factor for the third cycle (Tsang, 1983) The earlier used scheme which uses full penetration of the well in the aquifer during injection and production obtained a recovery factor of 0.-10 for the numerical simulation. By changing this scheme, a substantial improvement can be reached up to a recovery factor of 0.52.

A survey of international activities in aquifer thermal energy storage until 1981 was presented in Tsang (1931). He also mentions several topics which need further research. Among these are shape and location of the hydrodynamic and thermal front, thermal dispersion, diminishing the effect of a regional flow and thermal pollution.

If a full penetration of the well is used then in the aquifer the displacement of the water will be primarily horizontal. The thermal front between injected warm water and colder surrounding regions is then essentially vertical. However, such a thermal front is unstable due to the lower density of the warmer water. The thermal front will tilt. The hot water tends to rise and pushes the thermal front further from the injection well in the vicinity of the top of the aquifer. At the bottom, the hot water that has moved away will be replaced by colder water. This moves the thermal front towards the injection well. Especially aquifers with high

permeabi1i ty show this tilting. The effect during production is that besides hot water at the top also cold water at the bottom of the aquifer wil 1 be extracted.

The heat loss in the aquifer is roughly proportional to the area of the warm storage region and shall therefore be kept as compact as possible. However, a side aspect of this tilting is that this area increases. The rate of tilting of a more or less vertical thermal front in different aquifer types is studied theoretically by Hel Istrom et al (1979). They concluded that the rate of tilting is inversely proportional to the height of the aquifer. Furthermore, it is proportional to the square root of the product of the horizontal and vertical permeability. The rate of tilting due to buoyancy depends strongly on the injection temperature T, and the local temperature Tp, ■ It is increased five times, when Tg = 19. and T. = 99 . "C instead of 10. respectively 46."C. The diffuseness of the thermal front, however, diminishes the rate of tilting.

In connection with a field experiment carried out at Bonnaud (Jura) in France (Sauty et al, 1932b), a theoretical study on the thermal behaviour of a hot water storage system was made (Sauty et al, 1932a). The assumptions to simplify their numerical model of the description of this process are the following. Firstly, the aquifer is assumed to be homogeneous, of constant thickness and of an infinite lateral extent. Secondly, cap and bed rocks are assumed to be impermeable and of infinite vertical extent. Thirdly, the aquifer thickness is small enough so that vertical temperature gradients can be neglected in the aquifer. Fourthly, the aquifer and confining rocks are initially in thermal equilibrium with the geothermal gradient. Fifthly, the heat capacity and thermal conductivity are considered to be homogeneous and isotropic. The resulting mathematical model has two equations describing the heat balance: one for the confining layers (only diffusion) and one for the aquifer (radial diffusion and convection, and vertical exchange with the confining layers). They showed that if symmetric storage cycles (identical injection and production flow rates, and durations) are applied then the temperature distribution can be described with only four dimensionless parameters. The temperature distribution is also controlled by the number of cycles. The same derivation can be found in Doughty et al (1981). Sauty et al (1982a) also found that more cycles have to take place before the recovery factor remains constant if the heat loss to the confining layers is diminished. The production temperature is higher then. Depending on the

(15)

28

-r a t e of i n j e c t i o n and the heat loss to the c o n f i n i n g l a y e -r s , the f a l l of the p r o d u c t i o n temperature v e r s u s time i s f l a t i n the b e g i n n i n g and steep at the end of the p r o d u c t i o n p e r i o d ( 2T / t2 i 0) or r e v e r s e d ( 2T / t2 > 9> . The f i r s t behaviour occurs f o r the s i t u a t i o n when t h e r e i s a h i g h r a t e o-f i n j e c t i o n / e x t r a c t i o n and a v e r y low heat l o s s t o t h e c o n f i n i n g l a y e r s . The reason i s c l e a r , the area w i t h a u n i f o r m temperature equal to the i n j e c t i o n temperature has some e x t e n t t h e n . The second one occurs when t h e r e i s a low r a t e of i n j e c t i o n / e x t r a c t i o n and a high heat l o s s . The area i n the v i c i n i t y of the i n j e c t i o n w e l l , which has a temperature equal to the

i n j e c t i o n t e m p e r a t u r e , i s now v e r y s m a l l .

In S w i t z e r l a n d a major p r o j e c t named 5PE0S has been c a r r i e d out (Saugy et a l , 1983). The concept of u s i n g an a q u i f e r as s t o r a g e medium f o r s e n s i b l e heat i s r e a l i s e d i n a somewhat d i f f e r e n t manner. A v e r t i c a l w e l l - h o l e i s c r e a t e d i n the s u b s o i l . At two l e v e l s i n t h i s w e l l - h o l e , s i x h o r i z o n t a l d r a i n s are i n s e r t e d i n t o the s o i l . The aim of t h i s geometry i s t o reduce the e f f e c t of n a t u r a l convection by imposing a v e r t i c a l movement of the water due to a p o t e n t i a l d i f f e r e n c e between the upper and lower d r a i n s - The d i r e c t i o n of t h i s movement i s downward f o r the p e r i o d of i n j e c t i o n and upward f o r the p e r i o d of s u p p l y . However, f o r the p a r t i c u l a r s i t e chosen, a layer w i t h 1ow v e r t i c a l p e r m e a b i l i t y i s present between the upper and lower d r a i n s . The i n j e c t e d w a t e r , more or l e s s , s t a y s in the upper a q u i f e r and f l o w s o f f h o r i z o n t a l l y d u r i n g i n j e c t i o n . Natural convection due to buoyancy f o r c e s i s s i g n i f i c a n t because the v e r t i c a l d i f f e r e n c e i n p o t e n t i a l as i m p l i e d has a n e g l i g i b l e e f f e c t on the groundwater movement. The r e s u l t i n g heat l o s s to the atmosphere i s c o n s i d e r a b l e .

I n H o l l a n d a s t a r t has been made t o r e a l i z e an a q u i f e r thermal energy s t o r a g e i n c o r p o r a t e d i n a s o l a r h e a t i n g system of a l a r g e o f f i c e b u i l d i n g ( S n i j d e r s it De W i t , 1984) . In t h e i r study they conclude that a system w i t h h i g h performance c o l l e c t o r s o f f e r the p o s s i b i l i t y t o achieve much higher f o s s i l f u e l savings than the system w i t h the low performance c o l l e c t o r s . The f i r s t type of c o l l e c t o r s are a p p l i e d t o g e t h e r w i t h a s t o r a g e of high temperature l e v e l w h i l e the other one i s a p p l i e d w i t h a storage of low temperature l e v e l . To use the 1ow temperature l e v e l , a heat pump has t o be i n o p e r a t i o n . They conclude t h a t the low temperature system i s economical 1y more a t t r a c t i v e ,

CHAPTER 3 T^EOFSY ©F MATUKAL CONNECTION DM

P©FJ©US MEDIA

3.1 Introduction

Besides diffusion, convection of heat is one of the mechanisms of heat transport in a porous medium. Convection can be divided into forced convection and free (natural) convection. The latter is controlled by buoyancy forces. For a project of seasonal heat storage, these buoyancy forces result from the introduction of heat into the storage region of the water-saturated soil. Literature about the topic of natural convection is presented in this chapter.

The case that has been most extensively treated in literature is that of convective Darcian flow in a porous layer bounded by isothermal horizontal or vertical planes in the thermal steady state with the layer being infinite or finite in extent. The dimensionless Rayleigh number (Ra-number), modified for porous media, is introduced if there is a characteristic length H over which a temperature difference SÏ exists and if the parameters k, V, P are constant.

JT H K R3

=

Where K is the thermal diffusivity of the medium based upon the thermal capacity of the fluid and V and P are respectively the viscosity and the thermal expension coefficient of the fluid (water).

Besides this field also a porous layer with internal heat generation obtained some attention.

The main questions concerning natural convection are:

- the onset conditions for natural convection

(16)

- 30 ~ - the type o-f fluid -flow which occurs

3.2 Natural convection in a porous layer

The -first investigators (Horton & Rogers, 1945) dealt with the onset o-f convection in a horizontal porous layer con-fined between two in-finite sur-faces. They considered the solution -for steady free convection. The bottom face is impervious and isothermal. The fluid properties such as viscosity and thermal expansion coefficient were presumed to be constant. Lapwood (1948) studied the onset of natural convection when the upper face is impervious and isothermal and has a lower temperature than the bottom. He applied the linear stability theory. To perform this analysis, the dependent variables are considered to be the superpositi on of the solu ti on for the conduction only problem (stagnant flow and linearized temperature profile) and a disturbed part. Further it is assumed that the amplitude of disturbance (convection) is small. The resulting equation can then be

Table 3.1: Rar for several boundary conditions.

t o p K 9 1 8 n e l l 8 a l l 8 1 8 1 8 any bottom K e 1 8 1 8 1 8 1 1 1 1 1 8 1 1 8 any Rac A u t h o r ( s ) 4 K2 Lapwood 27.18 , , 27.18 1 N i e l d , \ Ribando 17.05 ) fk Torrance 12.8 N i e l d

'

2

,.

4- , , 8 . , , k=l impervious, rigid 0=1 constant heat flux k=0 constant pressure 0=0 constant temperature

linearized. Katto & Masuoka (1967) verified this value experimentally but showed that the thermal diffusivity, included in Ra, must be defined as the thermal conductivity of the porous medium divided by the thermal capacity of the fluid. Beck (1972), using the energy method, studied the effect of the positioning of vertical adiabatic boundaries on the onset of convection in an enclosed three dimensional porous medium. These lateral walls have little

31

-effect on the critical Ra-number Ra except in very narrow, tall boxes of which the width is much smaller than the height. Then Ra is larger than 4 it2- He also showed that the preferred cellular mode appears to be that the number of rolls and the direction of their axes are such that each roll has the closest approximation to a square cross section possible. The onset of free convection in a porous horizontal layer of infinite extent under other boundary conditions has also been discussed by Lapwood (1948) and Nield (1968) employing the linear stability theory. To study the onset of natural convection Ribando ic Torrance (1976) employed a finite difference method to

solve the governing equations. They found that in all cases Ra is less than 41(2 (Table 3.1).

Lapwood (1943) established that the critical Ra-number for impervious isothermal boundaries is also valid for the case of a free upper surface. He has shown that to a first order approximation the boundary conditions of a rigid surface and a free surface are equivalent in a saturated porous 1ayer .

To study natural convection with a Ra-number higher than the critical value, for which linear stability analysis is not valid, other techniques, such as perturbation methods, finite amplitude analysis, and numerical methods, have been used. Wooding ( 1957) used a per turbat ion expansion of the dependent variables in a power series in terms of a smal1 parameter proportional to

the Ra-number.

If the Ra-number is higher than its critical value, but not too large, a stable convectïve state exists. This convective f 1ow consists of adjacent hexagonal cells (or Benard cells) or of contra-rotating rolls which are two-dimensional. If the Ra-number increases and exceeds a particular value, the temperature distribution inside the porous layer begins to fluctuate continuously. This state corresponds with a relative increase of the heat transfer. H o m e & Ö'Sullivan (1974) interpreted this state as continuous creations and disappearances of convective cells, even in the thermal steady state.

In the previous investigations the assumption of constant viscosity and thermal expansion coefficient have been made. Wooding (1957) was the first who took the temperature dependency of the viscosity into consideration. He employed a perturbation method. The effect of a temperature dependent viscosity on the onset of cellular convection in a confined porous medium heated from below was studied by Kassoy & Zebib (1975). They employed the linear stability theory. The top and bottom temperature were respectively 25

(17)

32

-and 256'C. The variation in viscosity is more than a factor 3. The critical Ra number (with properties of the fluid evaluated at the upper surface temperature) was found to be substantially reduced. This is expected on physical grounds since the hot fluid near the bottom is considerably less viscous than elsewhere and hence becomes unstable sooner- than a relatively cooler fluid. Booker (1976) investigated the heat transport and f1ow

35 30 O Cfc 25 20 15 10 5

-_

"

-^Sfig\-\

?oaVf 200 . 1 , , . . , . , . v \ 25 \\so\

r

V7 5

-_

~

-\<

H (km) F i g u r e 3 . 1 : c r i t i c a l Ra f o r the onset of c o n v e c t i o n as a f u n c t i o n of the t h i c k n e s s of the porous layer f o r v a r i o u s values of T <"C/km) ( S t r a u s fe Schubert, 1977)

s t r u c t u r e of the temperature dependent v i s c o s i t y experimental 1y. In a d d i t i o n t o the temperature dependent v i s c o s i t y Straus & Schubert (1977) have taken i n t o c o n s i d e r a t i o n the temperature dependence of the thermal expansion c o e f f i c i e n t of water on the onset c o n d i t i o n s of f r e e convection i n t h i c k geothermal l a y e r s . They showed t h a t t h i s l a y e r i s much more u n s t a b l e than the l a y e r would be i f the water was t r e a t e d as a Boussinesq f l u i d having constant p r o p e r t i e s of s u r f a c e w a t e r . Also the c o n v e c t i o n i s c o n c e n t r a t e d i n the lower p a r t of the porous layer which i s h o t t e r . The major e f f e c t s are due t o the temperature dependency of the thermal e x p a n s i v i t y and the v i s c o s i t y . The c r i t i c a l Ra-number d e f i n e d on the thermal p r o p e r t i e s at the s u r f a c e i s reduced by 67. f o r a layer t h i c k n e s s as small as 18 meter. The bottom temperature i s 2D'C above the s u r f a c e temperature of 25*C. A c o n s t a n t p r o p e r t y c a l c u l a t i o n u s i n g mean v a l u e s f o r v i s c o s i t y and thermal e x p a n s i v i t y does not produce an accurate d e s c r i p t i o n of the f l o w f i e l d and temperature d i s t r i b u t i o n . The pressure dependence of the onset c o n d i t i o n s becomes not s i g n i f i c a n t u n t i l the temperature d i f f e r e n c e across the layer exceeds 256D'C and the temperature g r a d i e n t i s

everywhere g r e a t e r than 58rC/km.

33

-The i n v e s t i g a t i o n s mentioned thus f a r d e s c r i b e the onset of convection i n the v i c i n i t y of the c r i t i c a l Ra-number. I f the Ra-number i s increased beyond the c r i t i c a l v a l u e Ra then the amplitude o f c o n v e c t i o n grows and the t h e o r i e s based on a small amplitude of c o n v e c t i o n f a i l to d e s c r i b e the phenomena. Numerical methods are then used.

Elder <1967 a) showed by experiment and a numerical s o l u t i o n ( f i n i t e d i f f e r e n c e s ) t h a t f o r a system u n i f o r m l y heated from below, f o r which the f l o w i s c e l l u l a r , the Nu-number i s p r o p o r t i o n a l t o Ra-Ra f o r Rayleigh numbers (Ra> j u s t above the c r i t i c a l one ( Ra ) . The heat f l u x across the layer i s then p r o p o r t i o n a l to the square of the temperature d i f f e r e n c e across the layer but i s independent of the thermal c o n d u c t i v i t y of the medium or the depth of the l a y e r . For higher v a l u e s of Ra he found t h a t the Nu-number f o l l o w s the v i s c o u s f 1ow Ra-number w i t h power 1/4. T h i s occurs when the boundary layer t h i c k n e s s i s somewhat smaller than the c h a r a c t e r i s t i c l e n g t h of the s o l i d p a r t i c l e s .

For a v e r t i c a l s l o t f i l l e d w i t h a porous m a t e r i a l Nu i s a f u n c t i o n of Ra, the aspect r a t i o ( h e i g h t d i v i d e d by w i d t h ) and the boundary c o n d i t i o n s . The i n f l u e n c e of aspect r a t i o can be seen i n Bankval1 (1974) and Weber ( 1 9 7 5 ) . As c o u l d be aspected, i t i s found t h a t the heat t r a n s f e r decreases as the aspect r a t i o i n c r e a s e s . Th i s phenomenon i s e x p l a i n e d by the convec t i v e f 1ow from one v e r t i c a l s i d e to the other at the top and the bottom of the space. In a space w i t h a l a r g e h e i g h t the f l o w i n the end r e g i o n w i l l i n f l u e n c e a c o m p a r a t i v e l y small p a r t of the t o t a l space. Elder (1967 b) s t u d i e d the development of the f l o w n u m e r i c a l l y when a b l o b of hot f l u i d i s r e l e a s e d at the base of a hor i zon t a l porous s i a b or when a por t i on of the base of the porous s l a b i s suddenly h e a t e d . He also gave a d e t a i l e d comparison between the time-dependent s o l u t i o n s and the t i m e - l i k e development of the i t e r a t i v e s o l u t i o n s of the steady e q u a t i o n s . Q u a l i t a t i v e l y , t h e r e i s some aqreement between them.

During the l a s t decade several i n v e s t i g a t o r s solved a t h r e e - d i m e n s i o n a l model n u m e r i c a l l y . H o i s t & Aziz (1972 a&tb) were the f i r s t who gave the s o l u t i o n f o r t r a n s i e n t t h r e e - d i m e n s i o n a l f r e e convection i n a c o n f i n e d porous medium. I t was found t h a t w h i l e two-dimensional r o l l s t r a n s f e r more energy than t h r e e - d i m e n s i o n a l f l o w at low Ra-numbers (Zebib fe Kassoy, 1978), the r e v e r s e i s t r u e at h i g h e r Ra-numbers. Home (1979) n u m e r i c a l l y solved the governing equations which were f o r m u l a t e d w i t h the d e f i n i t i o n of a v e c t o r p o t e n t i a l . He found f o r the cubic r e g i o n t h a t the f l o w p a t t e r n at a p a r t i c u l a r v a l u e of Ra i s not unique but i s determined by the i n i t i a l

(18)

conditions. In some cases four alternatives exist: two- and three-dimensional, steady and unsteady flow.

3.2.1 Natural convection with internal heat generation

In a heat generating porous medium with impervious isothermal top and bottom where the bottom T. has a different temperature than the top T^ , the onset of convection is influenced by the rate of internal heat generation and the temperature difference between the bottom and the top. Gasser & Kazimi (1976) studied this interrelation which causes destabiIization CTu '■> T^ > or stabilization <TL < Tt ) of the fluid in

the porous layer. They considered the lower boundary to be impervious and the upper boundary to be a free surface.

t o p v. e 1 8 l l l e l i e e e 1 i e b o t t o m K e e e e ï ï e ï ï 1 e 1 ï 1 1 R ai , c 4 7 2 . 4 2 4 3 . 8 4 7 1 . 3 2 3 8 . 1 2 3 5 . 1 2 3 . 8 (51.9 A u t h o r i i k=l impervious, rigid 6=1 adiabatic

k=8 constant pressure 9=8 constant temperature i Kulacki & Ramchandani < 1975) ii Tveitereid (1977)

Kulacki & Ramchandani (1975) analytically determined the onset of convection in a horizontal porous layer with uniform heat generation by using the method of linear stability. As internal Ra-number Ra: can be defined:

g P k H M

Their results for several boundary conditions can be found in Table 3.2. Besides stability of the flows, hexagonal cells and two-dimen*"' onal rolls

35

-Tveitereid (1977) also investigated the onset of convection.

3.3 Some applications

Natural thermal convection in porous media occurs in a great number of actual situations. For some of these, convection exists and we are only able to observe it and try to understand and to control the phenomenon. On the one hand, this is the case for convection in geothermal areas, in oil reservoirs and in aquifers. The formation of geothermal reservoirs is believed to be intimately associated with the occurence of recent vulcanism or intense tectonic movements. As a result of these activities, magmatic intrusions may occur at little depths in the earth crust. Water is heated directly or indirectly by the intruded magma and is then driven buoyantly upwards to the top of the aquifer. If the aquifer is confined by caprock with a poor thermal conductivity which prevents the dissipation of heat then a large amount of hot fluid at relatively little depth could be available. Thus an ideal reservoir for geothermal energy extraction would consist of a hot heat source for a continuous supply of energy, a high permeable aquifer to permit the free convection of a large amount of groundwater and a poor heat conducting caprock for the storaqe of heat (White, 1973). A comprehensive review of the state of the art of the investigation of heat

transfer in geothermal systems is given by Cheng (1973).

On the other hand, in technical applications we are able to develop and use convection, e.g. by the storage of agricultural produce (Beukema, 1988), to control convection (thermal energy storaqe in an aquifer, the top layer of the soil or in a rockbed) or to reduce convection, e.q. for thermal insu1 ating purposes. It might be c1 ear that for our problem of using the upper layer of the water-saturated soil as a storage medium, we have to control and reduce the development of natural convection whenever it occurs

in order to reduce the heat loss.

From this investigation of literature in the field of natural convection in water-saturated porous media, it can be concluded that for an accurate description of this process by a mathematical model it is sufficient to incorporate density variations only in the terms where they create buoyancy. However, it is necessary to take the temperature dependency of the density and the viscosity of water fully into account. The mathematical model consists of a set of partial differential equations. These equations,

(19)

3é

-d e s c r i b i n g the - f l u i -d f l o w an-d t r a n s p o r t o-f h e a t , are n o n - l i n e a r an-d c o u p l e -d . To s o l v e these equations together w i t h the boundary and i n i t i a l c o n d i t i o n s , a numerical method has t o be employed.

37

-CHAPTER 4

THEORY AMD MATHEMATICAL MODEL

4.1 Introduction

In this chapter the describing equations for the problem of simultaneous fluid flow and heat transport through a porous medium are derived. Firstly, the concept of a representative elementary volume and the specific discharge which can be seen as an average velocity through a cross section of the medium will be given. After that, the proper equations describing continuity of mass, momentum and energy are presented for the water-saturated porous system. The main parameters which affect these equations will be discussed. The interaction between the storage reservoir and the energy supply and energy demand is represented by a source term in the energy equation and will be discussed separately. Finally, the boundary conditions imposed at the edge of the calculation domain and the initial condition are formulated. Then the problem is well posed.

4 . 1. 1 Representative elementary volume

Figure 4.1: The flow structure

Soil constituents exist in solid, liquid and gas states. The solid pha=f is represented by mineral and some organic particles with size variation -rom submicroscopic to visually discrete. The liquid phase is a solu' ,n of

(20)

38

-various salts in water. The gaseous phase consists mainly o-f air but with a high concentration o-f water vapour (Yong ii Harkentin, 1975). Solid

particles (clay, silt or sand) enclose void spaces which may be -filled with gas or liquid. Most o-f these voids (pores) are interconnected and form channels through which a -fluid may -flow; soil is a porous medium. The groundwater flows capriciously through the porous structure (figure 4.1). It is practically impossible to measure the velocity distribution in these flow channels. The flow equations applied at a point in the fluid together with knowledge of the complex structure of the flow channels with its irregular boundaries could give the flow distribution. This approach is termed microscopic and it may be clear that it will be a hard task to compute this flow distribution. A possibility to overcome this difficulty is to replace the actual ensemble of solid particles that make up the porous structure by a representative continuum. This can be done by setting up balance equations for a reference volume which is small compared to the dimensions of the porous medium, but large compared to the characteristic volume of the voids: Hubbert <1956), Whitaker (1966) and Slattery (1967). This approach is called macroscopic. Various properties such as pressure and temperature must be considered as averages over such a reference volume. This introduces what Bear (1972) terms a representative elementary volume (REV). Average values of dependent variables <PiT> and parameters (\, Pc ,n) can be assigned to mathematical points. These average values can be treated as continuous properties in space.

4.1.2 Specific discharge

Figure 4.2: The average linear velocity

Consider a cross section A of a part of a porous medium. H a volume flow

rate ♦ moves through this surface we can define the specific discharge

v . as:

This quantity has the dimension of a velocity and is also named Dsrcy- or superficial-velocity. The specific discharge can be related to an average linear velocity v. If the parameter n is defined as the porosity of the porous medium then n.A can be considered as the actual cross sectional area through which the fluid flow occurs. The relation between the specific discharge and the average linear velocity becomes:

This velocity v does not represent the true average velocity of the fluid elements in the pore spaces. This microscopic value is larger than v because the fluid elements must travel along irregular paths that are longer than the linearized path represented by v.

4.2 Flow equation

The basic flow description is represented by a system of equations stating conservation of mass and momentum. An Eulerian approach has been adopted for this description.

Similar to Bear (1972) and Freeze & Cherry (1979), the equations for the conservation of mass for a fluid flow through a porous medium on a differential volume in a general coordinate system can be written as:

3( nP ) ' f

To s i m p l i f y t h i s e q u a t i o n , the Boussinesq approximation ( B o u s s i n e s q , 1983) w i l l be a p p l i e d . T h i s approximation s i g n i f i e s t h a t d e n s i t y v a r i a t i o n s can

(21)

be neglected except where they create buoyancy. This is justified when:

(4.4)

The mass bal ance e q u a t i o n ( 4 . 3 ) -for an unconsol i d a t i n g porous medium

(n=constant) can now be r e f o r m u l a t e d a s :

Laminar -flow through porous media can be considered as c r e e p i n g f l o w . The c o n v e c t i v e con t r i b u t ion i n the equat ion f o r the conservat ion of momentum can then be n e g l e c t e d compared to the v i s c o u s f o r c e . As a r e s u l t , the d e s c r i p t i o n of t h i s balance e q u a t i o n f o r laminar f l u i d f l o w in i s o t r o p i c porous media y i e l d s D a r c y ' s law.

K

- ( grad p - P g ) <4.o) |t f

-T h i s law i s e m p i r i c a l l y determined by Henry Darcy ( 1 3 5 6 ) . Whitaker (l?oo) and S l a t t e r y (1967) d e r i v e d t h i s law from the N a v i e r - S t o k e s equation by-volume a v e r a g i n g .

S u b s t i t u t i o n of D a r c y ' s law equation (4.Ó) in the c o n t i n u i t y equation (4.5) r e s u l t s i n the f l o w e q u a t i o n :

d i v < - < grad p - p g > ) =

\i 4

-We mainly considered the deviation p' of the hydrostatic pressure during our computations. This hydrostatic pressure is based on a reference density Pa ■ The flow equation now yields:

k

div < - [ grad p' - < p - p > g ] ) = 8. (4.8) H t 8

-The p r e s s u r e d i s t r i b u t i o n and the groundwater f l o w p a t t e r n are coupled w i t h the temperature d i s t r i b u t i o n through the d e n s i t y P<T) and v i s c o s i t y t»(T) .

41

-At each time t h e r e i s a c e r t a i n temperature d i s t r i b u t i o n i n the s o i l . T h i s d i s t r i b u t i o n determines the development of n a t u r a l c o n v e c t i o n . The temperature changes w i t h time due t o , among o t h e r s , the c o n v e c t i v e and d i f f u s i v e heat f l o w s . The buoyance f l o w w i l l g r a d u a l l y change w i t h time due to the changed t e m p e r a t u r e , d e n s i t y and v i s c o s i t y d i s t r i b u t i o n s throughout the s o i 1 .

Given a temperature d i s t r i b u t i o n , the d i s t r i b u t i o n of the d e n s i t y and v i s c o s i t y can be c a l c u l a t e d . The f l o w equation ( 4 . 7 ) can now be r e s o l v e d to g i v e the pressure d i s t r i b u t i o n . A f t e r w a r d s , D a r c y ' s law (4.Ó) can be used t o determine the d i s t r i b u t i o n o f the s p e c i f i c d i s c h a r g e .

4 . 2 . 1 L i m i t a t i o n s of D a r c y ' s law

The d e c i s i o n t o use D a r c y ' s law as the e q u a t i o n f o r the c o n s e r v a t i o n of momentum depends upon the dimension of the REV (1) and the c h a r a c t e r i s t i c l e n g t h d . For l>>d , the porous medium may be considered as a continuum and D a r c y ' s law may be used. However, when 1 and d becomes of the same magnitude, the assumptions of WhitaKer (19óé) are no longer s a t i s f i e d and the v i s c o u s term may not be approximated by the Darcy-term i n the volume averaged N a v i e r - S t o k e s e q u a t i o n . An upper boundary f o r D a r c y ' s law i s where the Reynolds number ( v . d / v) becomes around 1. T h i s Reynolds number r e p r e s e n t s the r a t i o between the c o n v e c t i v e and v i s c o u s f o r c e s . Hubbert ( 1 9 5 4 ) , among o t h e r s , r e c o g n i z e d t h a t D a r c y ' s law i s no longer v a l i d a t h i g h r a t e s of f l o w . At t h i s p o i n t the f l o w i s s a i d t o be n o n - l i n e a r because the

Reynolds number. Re

. io~z tcr' TO JO' IO2 X)3 CD L lineor_ _Jk>^r\ _l u r b u l e n l jJj | laminar Uammar]

Figure 4.3: The specific discharge as a function of Reynolds number

convective force may no longer be neglected. The characteristic length in Re is not a unique value. In literature several parameters such as mean

(22)

42

-pore d i a m e t e r , mean p a r t i c l e diameter and the square r o o t of the p e r m e a b i l i t y are taken f o r t h i s l e n g t h . Bear (1972) summarizes the experimental evidence w i t h the remark t h a t D a r c y ' s law i s v a l i d as long as the Reynolds number based upon the average g r a i n diameter does not exceed a v a l u e between 1 and 16 ( f i g u r e 4 , 3 ) . 4 . 2 . 2 P e r m e a b i l i t y

-Glacial

tilt--Silt,

loess--Silty

sand-Clean sand—

Gravel-1Ö

20

10'

8

10

,e

tO 10'*

2

10-'° 10

s

k(m

2

)

Figure 4,4l The permeability of several types of soils

The intrinsic permeability V, (m2) can be defined as the ability of the porous medium to transmit water. It is only a function of the porous matrix. The matrix properties which are relevant to the intrinsic permeability are e.g. the porosi ty, the grain size distribution, the shape and texture of the grains and their specific surface. The permeability for several types of soils can be found in figure A.A.

For anisotropic porous media Bear (1972) developed a similar expression for v , as equation (A.i) in terms of a permeabi1lty tensor k. In the systems

we wish to model little is known regarding the tensorial properties of k. Further we think that large scale inhomogeneities will dominate the convective behaviour for the cases we are interested in. This is the reason that k will be treated as a scalar which can vary spatially. However, the model that is given above can describe anisotropy when the principle directions of anisotropy are parallel to the coordinate directions. The diagonal entries of the permeability tensor are then non-zero and the off-diagonal entries are zero. The resulting equation (4.7) is then the flow equation for laminar fluid flow in anisotropic, heterogeneous,

unconsolidating porous media.

A.2.3 Streamfunction formulation

For two-dimensional f l o w problems i t i s o f t e n convenient t o change over t o the streamfuncton f o r m u l a t i o n i n s t e a d of a d o p t i n g the pressure as dependent v a r i a b l e . The streamf u n c t i o n f can be de-fined a s :

SI 'óz

In this formulation, u and w are respectively the x- and z-component of the specific discharge. This definition satisfies the continuity equation. To find the streamfunction formulation of the flow equation, the pressure has to be eliminated from Darcy's law by cross differentiation of the two components of the specific discharge. The result yields,

ii a* k 3z 11 3¥ k 3x

_

3

^

3x

when the positive z-direction is pointed downwards.

A.2.A The density and viscosity of water

t_

temperature (°Ci

Figure 4.5: The density and viscosity of water