Speleogenesis and Evolution of Karst Aquifers

The Virtual Scientific Journal

ISSN 1814-294X
www.speleogenesis.info

Underground drainage systems and geothermal flux

Giovanni Badino
Dip. Fisica Generale – Università di Torino, Associazione La Venta
E-mail: badino@to.infn.it
Republished from : Acta Carsologica 2005, 34 (2), 277-316

Abstract
The paper presents an analysis of the interaction between the geothermal flux and the water or air- deep drainage
networks. The problem of geothermal power intercepted by deep structures and, in general, the temperature field
calculations, is converted to classical thermo-engineering problems in terms of shape factors. It is shown that the fluid flow
in a conduit perturbs the whole deep rock temperature field until the geothermal flux of a large area is focalised onto the
conduit. It is shown that either small water masses flowing into a mountain are able to perturb the rock temperature up to
the surface, on sizes that do not depend on water mass dimension, but on its depth, and then on enormous volumes. The
introduction of the “geothermal cross section” of an underground drainage structure allows us to improve the classical
formula of minimum provenance depth of geothermal water. Enlarging factors are applied to the classical estimation in
dependence of the ratio between the actual average discharge and the critical discharge Qc, which depends on the conduit
geothermal cross section. The geothermal “umbra cones” created in the overlying rock by deep underground structures are
described.
It is shown that the geothermal flux can play a significant role in the underground drainage phenomenology.

Keywords: geothermal flux, heat transfer in karst massifs

1. The geothermal flux and the energy The flux is some 60 kW per square kilometre,
contents of rocks which corresponds to a total release of 3×1013
W on the whole planet. For comparison the
An introduction to the geothermal energy flux energy flux received from the Sun is 1.7×1017
W, therefore the geothermal flux is around 5000
It is widely known that the rocks below us times smaller than the main energy source for
have temperatures that increase with depth. The Earth. It cannot play a role in the free
reason is that the internal part of Earth is hot atmosphere phenomenology, but we are going
and the surface cold; there are then two “heat to see that in the case of underground
sources” (but in this work I have adopted the atmospheres it can and it does.
suggestion of Bohren (Bohren, 1998), avoiding
use of the word “heat”), and the thermal energy At first, this appears not to be true in the case
flows between them with the rules given by the of cave atmospheres, that are really quite cold
usual conduction equations. (essentially at the yearly average external
temperature, -from hereafter Tave ; Badino,
Table 1, adapted from (Lee, 1966) gives 2004) while the atmosphere of mines can be
typical values, widely variable, of geothermal very hot (Badino and Forti, 2005). Actually the
flux, estimated by measures of deep first aim of this work was to show that the
underground temperature gradients. geothermal flux could not play any role in the
The world average (Verhoogen, 1956) is deep karst microclimates and genesis, because it
is shielded by deep drainage conduits: This is
Fgt = 0.06 W m −2 exactly the contrary of what we are going to
show...
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.2

TABLE 1
Typical values of geothermal flux
Regions Mean geothermal flux Fgt Number of measures
[mW m-2] or [kW km-2]
Continents
Africa 36-61 15
America 25-150 85
Australia 35-160 65
Asia 22-150 60
Europe 26-140 60
Oceans
Atlantic 3.4-250 250
Indian 5.9-220 250
Pacific 2.9-340 600
Arctic 33-62 20
Geothermal regions
Larderello (I) 450
Oora Kei Korake (NZ) 4200
Matsukawa (J) 630

Mines, caves and tunnels The Gran Sasso rock is essentially limestone,
The first point to discuss is the difference and drillings have shown a cave layer 550 m
between the temperature of artificial and natural above the tunnel altitude, and a general water
cavities in order to clarify the common table extended up to that level. The infiltration
confusion between “thermal flux” and waters at the surface are essentially at Tave, and
“temperature”. cross the whole mountain in nearly adiabatic
conditions, which means that they are only very
It is useful to compare two Underground slowly heated along the fall. Therefore in the Mt
Neutrino Observatories, in Mont Blanc Blanc’s depth there are essentially “mine”
(between France and Italy) and in Gran Sasso waters in thermal equilibrium with hot, deep
(in central Italy), which are assembled in halls rocks, whilst the deep Gran Sasso waters are
in motorway tunnels. The depth of the first, essentially meteoric waters, in equilibrium with
dismissed in 2001, was some 1800 m below the the atmosphere.
surface, at an altitude of 1300 m asl, whilst the
depth of the second is around 1050 m at an It is useful to discuss a little more the internal
altitude of 1000 m asl. Their temperatures are water heating in karst. The reason for water
nevertheless completely different, in the first it temperature increase during underground fall is
is around 32 °C, in the second at 6 °C, the two the gravitational energy which is converted in
unsuitable for working. Therefore it is necessary thermal (a process that gives a water adiabatic
to act in the opposite sense, in the first to cool, lapse rate –2.34 °C km-1) and, in the case of
in the second to warm the experimental halls. flowing in vadose conditions, also to thermal
exchanges with moist air, always characterised
The reason for the two different temperatures by a different adiabatic lapse rate, around -5 °C
is the different rock permeability. km-1. The actual caves’ lapse rates are between
The Mont Blanc rock is mainly granite, the these two extremes (Badino, 2000; Luetscher
waters met by the tunnel were essentially fossil and Jeannin, 2004), a fact that has huge
waters, the meteoric water circulation being consequences on the caves energetic balances,
quite epidermal (up to depth of 100-200 m which nevertheless are outside this work aim;
below the surface), with some exceptions we return briefly to it in the next chapter, but a
localised along large major rock discontinuities. discussion can be found in (Badino, 2005). We
concentrate here on the fact that the infiltrations
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.3

create a thermal connection between the

atmosphere and the deep karst.
In short the thermal contact between deep
waters and atmosphere in the Mont Blanc case
is due to the recent tunnel, in the Gran Sasso to
the natural ancient conduits that have permitted
a deep water flow that has been able to shield
the geothermal flux from below and to cool the
whole mountain above the water table.
More in general we can compare the
temperature in the bulk of large mountains Fig. 1. Depth profile and rock temperatures in the St.
looking at the temperatures encountered during Gotthard tunnel.
the tunnels construction. The world’s deepest
tunnels are essentially in the Alps and it is
possible to discuss their internal temperatures.
The figures, adapted from (Szechy, 1973) show
the rock profile above the tunnels and the
corresponding local temperature.
The first figure (Fig. 1) shows the situation of
the St. Gotthard tunnel, in Switzerland, mainly
in granite and gneiss. It is possible to see that
the temperature dependence with depth is quite
regular.
Next figure (Fig. 2) shows the situation in the Fig. 2. Depth profile and rock temperatures in the
Simplon tunnel (between Italy and Switzerland), Simplon tunnel.
a geological structure in gneiss and, roughly in
the Italian part, limestones. It shows low
temperature anomalies in the sedimentary part
around PK 15 (Luetscher and Jeannin, 2004).
The Mt Blanc situation is quite different (Fig.
3). Its rock is mainly protogine, but there is an
important fault that allows very deep glacial
water circulation, which has lowered the rock
temperature near PK 8, where extreme
excavation problems where encountered
(Guichonnet, 1967). The general behaviour is
regular, but the whole mountain has been cooled
a little by the fault. Fig. 3. Depth profile and rock temperatures in the
Mont Blanc tunnel.
The next figure (Fig. 4) shows a completely
different situation, the Gran Sasso (Catalano,
1993). Not only the temperature does not
increase with depth, but also it shows a
tendency to decrease with it, because as deep as
it is, as colder are inflowing waters. We have
seen above that they meet a very small warming
crossing the mountain (Badino, 1995).
Other more complex phenomena can occur in
determining the rock temperature. For example
Szechy cites the case of the Great Appennine
tunnel (Italy), mainly in limestone and clay,
where a local temperature of 60 °C, 500 metres Fig. 4. Depth profile and rock temperatures in the
Gran Sasso tunnel.
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.4

below the surface, has been encountered. The cannot be extremely dangerous. Otherwise, if
Author says that “discrepancies are due to during excavation appears a tendency to a
intensive methane infiltration through the temperature reduction, it is a sure signal of an
Eocene clay shale”, which looks quite strange. approaching water stream in direct contact with
Nevertheless it is quite common to meet high the surface, that at these depths has obviously
temperatures connected with local thermal water enormous pressure, which is able to create
infiltrations, often connected with hyperkarstic extremely dangerous situations, also because it
phenomena (Badino and Forti, 2005). Here we is surely associated to dramatic rock
deal only with a “standard” situation to analyse discontinuities. This is the reason why a lot of
the geothermal impact on our usual karst. work is made in the field of rock temperature
estimations in deep tunnels (Koenigsberg, 1906;
Convective mountains Goy et al., 1996; Badino, unpubl. 2005).

The graphics (Fig. 5) show in short the above Underground high temperatures are
discussed data (temperature vs. depth) where connected with good thermal insulations, which
each type of data point indicator describes a means that a tunnel, or a mine, can be excavated
different tunnel. It is possible to see two across rocks that are very hot (high temperature)
completely different behaviours, the “hot” because, i) they have almost no contact with the
mountains without internal water fluxes, with surface and ii) they have then acquired
positive temperature gradient, and the “karstic” equilibrium with the geothermal flux. Then the
mountains with slowly negative or zero relative rock insulation has allowed to a so
gradients. small thermal flux to heat up to high
temperatures enormous quantities of matter.
As we have seen above, it is possible to
consider mountains with caves as good thermal
conductors and then in general (let us forget for
a while the geothermalism) they are in thermal
equilibrium with the external atmosphere.
As larger are the caves, as smaller it is the
impedance for water crossing the whole
mountain, then deep water circulation is
possible and the whole structure is crossed,
which is a big difference from a mountain
without caves where only water circulation is
Fig. 5. Rock temperature versus depth in the large
quite epidermal.
alpine tunnels. Another consequence to be pointed out is that
if the mountain is highly permeable to water
fluxes the thermal energy transfers inside it are
In fact these mountains are in thermal absolutely dominated by fluid motion, then the
contact, and in equilibrium, with the local pure conductive terms into the rock (“heat”
climate, that they follow with a delay that
diffusion) become negligible.
depends on the mountain depth (Badino, 2004).
It is possible to add that low-impedance water
This thermal contact surface-underground (or air) transfers inside a mountain can be
due to water vein, is the base for the traditional considered a special case of thermal convective
and fundamental “remote sensing” during tunnel movements, and then the karstic mountains are
excavations. Continuous rock temperature examples of convective thermal contact with the
measurements are performed during work atmosphere, whilst the impermeable rocks have
(Guichonnet, 1967); a regular temperature only a conductive contact with it, with a very
increase with the surface distance is a signal of
poor efficiency.
compact rock. Occasional water veins are in
thermal equilibrium with geothermal flux,
The hot mountains energy contents
which is possible only if these water reservoirs
are relatively small and without a hydro- Let us estimate the energy contents of a hot
geological connection on long ranges, and then rock prism, with a surface A and an altitude of
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.5

H. Its temperature at the beginning is T0=Tave, power, but it has worked from a very long time,
local average temperature. the rock is a very efficient thermal insulator and
the final result is an efficient energy storage.
The geothermal flux will heat the rock until
the temperature gradient in it becomes able to
evacuate the geothermal power on the surface. Rock heating time
This equilibrium condition (“stationarity”, i.e. To estimate correctly the time needed for
no time dependence) is obtained when the heating, it is necessary to take into account that
gradient is the rock prism is in contact with others all
around; the calculation would have to include
∂T
Fgt = K R these in the estimation, reconstructing the whole
∂z temperature field and its dependence on time.
The problem is complex but we can do a crude
Where we have used partial derivative estimation of the prism heating time scale
because in general the temperature is function considering it as thermally isolated from its
also of time. This is exactly the condition that surroundings (which is equivalent to assume a
allows the measure of geothermal flux from flat surface and uniform, uni-dimensional
temperature gradient data, assuming a steady- heating), then
state condition.
In this case we can assume a rock thermal Q  Fgt Z 2  1
∆t heat = =  A C R ρ R  =
conductivity KR=2.5 Wm-1 K-1, typical for AFgt  K R 2  AFgt
granite, whilst the limestone conductivity is
some 10% less. Then the temperature gradient C ρ Z2
at the equilibrium is
= R R
KR 2
∂T Fgt 0.06
= = = 0.024 K m −1 The last term can be rewritten in terms of
∂z K R 2 .5 thermal diffusivity coefficient defined as
Which is 24 °C per kilometre. It is easy to KR
calculate the energy necessary to heat at this aR =
constant temperature gradient a rock prism of CR ρ R
surface A from depth Z up to the surface; its
total available energy is In the case of rock

Fgt Z 2 aR ≈1.2 × 10−6 m 2s −1

∫
Z
Q= AC R ρ R ∆T dz = AC R ρ R
0 KR 2 And then
Where CR is the rock specific thermal Z2
capacity (800 J kg-1K-1) and ρR its density ∆theat =
(roughly 2600 kg m-3). For example, assuming 2 aR
Z=1 km, with typical values the total energy
contents per kilometre square of surface is Typical width of limestone mountains are
around 1000 m, then the heating time scale is
Fgt Z 2 around 104 years, not so much for geological
Q = AC R ρ R = time scale, but longer than the typical global
KR 2
climate fluctuation (Badino, 2004).
1000 2
= 10 6 × 800 × 2600 × 0.024 × = 2.5 × 1016 J
2 Penetration lengths of temperature fluctuations

For comparison, a 20 kTon atomic bomb Let us recall the classical thermal fields
yields 1014 J. solution of a homogeneous thermal conductor to
a sinusoidal and to a sudden (step) temperature
It is important to note that this energy, that a change.
power plant of 1 GW produces in two years, has
been released to our prism rock by the small In the first case a thermal wave propagates
geothermal heater, only 60 kW. It is a small inside, fading exponentially (Badino, 2004;
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.6

Lismonde, 2002). If tsin is the period of temperature increase at depth x at time t, by

sinusoidal wave we have that the “penetration equation (Isachenko, 1969)
length” is
  x 
aR T (x,t ) = ∆T 1- erf  
l p2 = tsin   2 aR t 
π
If the temperature fluctuation has a step Where erf(u) is the Error Function. The figure
shape, the propagation length is more difficult (Fig. 6) shows the results at different times. A
to define. The penetration of thermal shock of discussion about the properties of this solution
amplitude ∆T is described in terms of rock can be found in (Lismonde, 2002).

Fig. 6. Diffusion in rock of

one step temperature
increase on its surface.

The last equation suggests that a step thermal temperature of the surface, but when its
wave is able to reach a depth x in a time t or, temperature has attained the “stationary
vice versa, after a time t the thermal disturbance temperature field” geothermal gradient seen in
has reached the depth x, and the relation the previous chapter. Nevertheless it makes no
between the two quantities is given by the sense to try to perform exact calculations, which
argument of Error Function. More precisely it is in any case work with unrealistic forms of
possible to show (for example, classically mountain.
(Laidler and Meiser, 1995) that We can conclude that the heating and the
cooling time-scale up to the equilibrium
x 2 ≈ 2aRt configuration are almost the same, and they do
not depend on the temperature drop, then at
Then the penetration in rock up to depth Z of
depth Z in metres we have
a cool wave requires a time-scale ∆tcool
Z2
Z2 ∆t eq = ≈ ∆t cool ≈ ∆t heat [1.1]
∆t cool = 2a R
2 aR
In the case of compact rock we can assume
The estimation is heavily approximated, that the equilibration time scale in years is given
because we are not looking for the complete by
cooling of the mountain, but for the equilibrium
temperature field formation inside it, what is ∆t eq = 0.01 Z 2 [y]
attained not when the whole rock is at the same
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.7

This estimation is very important for our complex, but it is better to postpone a detailed
discussion. In the next chapters we are going to discussion to a future work.
consider water and geothermal fluxes that have We can spend some words about the exact
attained stationary conditions. It is then obvious value of T0 and its relation with Tave. Really Tave
that if the water flow has begun from a time that depends on the altitude, and on average
is much smaller than that of equilibration time, decreases of 6-6.5°C per kilometre, as described
it is not possible to assume that the system has with the International Standard Atmosphere.
attained a general equilibrium configuration, Also at a first approach it looks better to think
what in fact happens in the artificial that the rock assumes not the average yearly
excavations. This is the reason why the mines temperature of the atmosphere Tave, but the
are hot. average temperature of waters at the infiltration
point, that is quite lower because the rain waters
in alpine karst are generally associated with
2. The basic approach to the problem colder periods (but in tropics with warmer
periods); really, many other corrections are
The problem of infiltration temperature necessary to estimate the local rock
temperatures (Badino, 2004).
The classical Desio’s formula (Celico, 1986) Another corrective term, already cited above,
that gives the minimum depth attained by appears during the underground flow because
geothermal water is based on the assumption the temperature increase of underground waters
that underground there exists a first layer along their travel is different from outside,
(heterotermic) in thermal contact with the where the waters follow the ISA mean lapse rate
atmosphere with essentially the temperature of (6 °C per kilometre of fall). Underground, in
local inflowing fluids, which have a temperature adiabatic conditions the theoretical value of
T0 quite precisely (but in general, a little lower, water temperature increase is 2.34°C per
(Badino, 2004)) equal to the local yearly kilometre of fall, but the experimental values
temperature Tave. Below this first layer the rock into the caves, where energy exchanges with the
temperature increases in conductive regime air are possible, are around -3 and -4°C km-1.
(homotermic region). Really, we can assume that the water
The heterotermic layer is also called “active temperature that has infiltrated at altitude H
layer” (US Bureau of Mines, 1996), mainly [km], when it arrives at level 0 is some (3-4)×H
because the seasonal variations can create icing. [°C] hotter, then sensibly colder than the
It is necessary to note, nevertheless, that corresponding infiltrating waters at that altitude.
some problems do exist in the definition of In detailed calculations these effects, that
homotermic layer (Schoeller, 1962; Celico, create a difference between the actual cave
1986), because it is considered the layer where temperature T0 and the local yearly temperature
temperature does not depend on atmospheric at its same altitude Tave, have to be taken into
temperature variations, not the layer where account, but in our discussion they are
temperature is equal to the local average of the completely negligible.
atmosphere. It is then defined on the basis of its
temperature stability, not on the basis of its The energy release to groundwater
thermal contact, if with atmosphere or with deep
rocks. Really in rock with deep aquifers we There is an obvious approach to consider the
meet a heterotermic layer with seasonal role of geothermal flux.
fluctuations (some score of metres), a layer In the upper Earth surface layers the
above the aquifer (included) at T0 (Luetscher, geothermal energy is essentially intercepted by
2004; Badino, 2005), a thin layer of thermal water that releases it to the atmosphere when it
contact deep rocks-aquifer, where a relatively goes out from springs.
sudden temperature increase, probably The energy that comes onto a large surface A
dependent of aquifer permeability, is possible is obviously FgtA, and it is very regular in time.
(Goy, 1996; Badino, unpubl. 2005). We would Let us suppose that it is absorbed by a mass M:
call it the “geothermal exchange layer”. Below Its temperature increase rate is then given by the
this region, a regular temperature increase in condition
deep rocks is found. The scenario is then more
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.8

Fgt Adt = Cw M dT Let us estimate the water flux out-flowing

from a region of surface A, which thermally
Where Cw is the water thermal capacity, that interact underground. If the precipitation is P (in
is kg m-2s-1), the infiltration is P minus the part Pout
“lost” outside due to evaporation. This part
dT F A depends on climate, surface type and so on; in
= gt temperate regions ranges between 30 and 40%
dt Cw M of total, but in deserts can rise to 90% (Celico,
Then if we know the shielding mass M we 1986). With this assumption the enthalpy
can calculate its temperature variation with extracted from the system is (P-Pout)∆TgtA,
time, but how large is M? Reasonably it is the where ∆Tgt is the water temperature increase
mass of groundwater, but it changes seasonally, during deep flow (Fig. 7). The condition of
depends on rock characteristics and so on. Still stationarity implies that on average the
worse, the previous equation says that the temperature cannot change in time and then
temperature continues to increase because the Fgt A = (P − Pout ) C w ∆Tgt A
mass is stable and the energy flux continues.
The equation gives us a temperature-changing That solves our problem. We can change
rate, not a temperature variation. units, calling P* the infiltration in [mm a-1], to
It is easy to avoid these problems. Really the obtain
mass M can change, but its water is almost
continuously flowing, which means that some
0.06 500
water enters cold in the mass M, and on average ∆Tgt = = [°C] [2.1]
the same quantity flows away from some other 4.2 × 10 (P − Pout ) P *
3

side, warmer. We will impose a condition of

steady state and calculate the temperature
variation of the flowing water, independently on
M and, above all, on time.
The geothermal energy flux is very regular
but the groundwater flow in general it is not.
We then may expect that the temperature
change fluctuates, also if we consider the over-
simplified system of a single water mass M and
we neglect other problems like the drainage
network structure and mixtures between
different branches that depend on the water flow
rate.
So, the temperature changes fluctuate, but we
are going to estimate the average value of
temperature increase due to geothermal flux.

The groundwater temperature increase

It is easy to estimate the flowing water
Fig. 7. Interception of geothermal energy flux by a
temperature increase, assuming the powerful flat aquifer. In stationary approximation the water is
and reasonable hypothesis that on average the heated and the upper rocks are completely shielded.
system is stationary. This means that the
thermal energy inflowing in M from Earth is, on
average, going out as water flux enthalpy The upper parts of drainage systems (for
increase. In this way we have not to consider the instance the caves, highly permeable) are almost
mass M, that we cannot know, but only the exactly at the external average local
outgoing flux from it, because our steady state temperature, therefore the water is in general
assumption states also that the temperature of M warmed of ∆Tgt between the lowest cave parts
does not change with time (on average...). and the springs, that is along the flow in the
phreatic systems.
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.9

In alpine karsts P* is some 1000 mm a-1 and assumptions true for similar “discrete” systems?
then the water average temperature increase due The general answer is that no, they are not
to geothermal energy is some 0.5°C. generally correct.
It is a small term, very difficult to measure, The water flows along definite branches, that
and really it has never been measured. But it is do not cover a large surface and with a general
not always so small, a temperature increase of shape that is far from regular. The complete
5°C between the caves and the springs has been shielding assumption is not reasonable for
measured in Cuatro Ciénegas, a figure that karstic drainage.
corresponds to an infiltration (P*-Pout)=100 mm
a-1 in this desert region of Coahuila, Mexico
(AA. VV., 2004). 3. The underground temperature field
with a drainage network
Therefore it would theoretically be possible
to estimate the average infiltrating flow
measuring this temperature increase, but it is a The problem
measure quite difficult to do with some The scenario described above (a regular, flat,
accuracy, because it is the difference of two diffuse water table) can sometimes be correct
uncertain parameters, and other processes surely but in general it gives completely wrong results
occur along the water rock crossing. in the internal rock temperature field
estimations.
The plane watertable
To study a more real model it is necessary to
We have used these calculations in a previous estimate the energy interception made by a
work (Badino, 1995) to explain why the karstic system (a thermodynamical sink) that is buried
mountains are so cold, therefore excluding a in a semi-infinite medium where a thermal
significant role of geothermal flux in energy flux is flowing from infinite.
underground climate definition, because the
Let us consider the problem details.
energy flow from Earth depths is easily
shielded. This very natural approach to refuse, We have a semi-infinite rock volume in
in general, any role to the geothermal flux, it which a thermal flux Fgt is coming from below.
was not original, it is a quite traditional point of It creates a temperature gradient given by:
view.
∂ T Fgt
Bögli (Bogli, 1980) estimates reasonably 0.2 =
°C “to prevent the karstified zone above from ∂ z Kr
being geo-thermally heated up”, that is a very Where we have assumed the depth z as
interesting idea that we shall meet again. In a positive downwards. The temperature field at
very interesting and complete paper Mathey depth H below the surface (or, better, below the
(Mathey, 1974) estimates a maximum of 0.55 heterotermic layer) it is therefore given by
°C. More recently Jeannin et al. (1997) estimate
the specific discharge of a karstic spring (the ∂ T  F 
equivalent of infiltration) between 30 and 3000 T (H ) = T0 +  H = T0 +  gt H
mm a-1. The first figure seems too low (almost  ∂z   Kr 
three times dryer than the North Mexico
We have previously seen that T0 is essentially
deserts), but in this paper it is used an energy
the local Tave. It is very natural to suppose that
flux too low by a factor thousand, to finish to
the geothermal energy intercepted by a deep
say that the temperature increase, that has
structure is that given in previous chapters, that
physical dimension [°C], is “less than 0.1 °C/a”,
where really the 30 mm a-1 case would be is (geothermal flux)×(structure area), therefore
warmed by some 15 °C. W = AFgt [3.1]
Let us return to this estimation. We have
obtained the average water temperature increase This means that the energy intercepted by a
and we could now begin to calculate the “cave” is proportional to its surface in the
fluctuations dependence on discharge and so on. direction of energy flux. It is a very natural
But we are dealing with conduits, not with plane assumption, but it is false.
watertables. Are these calculations and
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.10

The fact that, up to now, has not been the stationarity assumption. Whatever initial
considered is this: If we bury a system able to temperature field condition will converge to
intercept and to evacuate elsewhere the energy, asymptotic values which are solution of Laplace
the whole temperature field in the rock is equation, but this convergence requires time.
altered and then the flux itself changes. During this time, which is of the order of
equilibration time scale introduced above [Eq.
The geothermal field with a cave 1.1], the difference between the real field (in
transient phase) and the asympotic one
The problem of temperature field calculation (stationary phase) can be important. If the
in this configuration has to be solved with these equilibration time is comparable with the typical
boundary conditions: changes of boundary conditions (global climatic
1) The geothermal flux from infinite is changes, infiltration of hot waters or so on), the
constant, Fgt; system can never be considered in a true
2) The temperature on the surface is T0 stationary phase, and the equation given by this
constant everywhere; assumption has to be considered heavily
approximated. This limit of stationary
3) The cave temperature is T0, the same as on
assumption gives strong uncertainty in the
the surface (we have sufficiently discussed the
temperature fields estimations for the new deep
limits of this assumption);
alpine tunnels, which have very long
4) At the infinite the temperature field is not equilibration times (Badino, unpubl. 2005), and
disturbed by the cave existence. affects also our next considerations.
These conditions imply a flat external
topography and assume that the infiltrated water
heating from surface to the cave is negligible.
We said that if the cave has an area A, it is
natural to assume that the absorbed geothermal
energy is Fgt×A [Eq 3.1]. Is this correct? Let us
firstly discuss it qualitatively, drawing the
isothermal surfaces in the rock. The figure (Fig.
8) shows a reasonable situation that respects the
boundary conditions. It is possible to see two
things:
1) The isothermal surfaces have a tendency to
converge, then to be focussed, onto the cave;
2) They are “compressed” around the cave.
If we remember that the thermal flux flows
along the maximum T variation (i.e. along the
grad (T), which means perpendicularly to the
isothermal surfaces) and that its value is Fig. 8. Qualitative view of the stationary thermal
field due to the interaction between the geothermal
proportional to the gradient of T, we have that
flux and a conduit with strong drainage (system S).
the two features are equivalent to say that: i) the The isothermal surfaces are affected in a wide
cave focus on itself the geothermal flux and, ii) region, and the cave temperature is T0.
in the rock surrounding the cave the geothermal
flux (and the geothermal gradient) is much more
intense than the natural one...
The general solution
So, the assumption that gave us the [Eq. 3.1]
is surely wrong. Therefore the problem of calculation of
stationary temperature fields is very difficult to
But is it possible to calculate the correct
solve just with the easy boundary conditions
value? Quite surprisingly it seems that nobody
given in the previous chapter. What to expect,
has studied this important problem.
then, when we will have to assume finite energy
Before we look for the solution, we have to transfer rate inside the cave or situations in
make some note about the enormous weight of which the temperature of the cave itself is
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.11

determined by the geothermal flux? And what where H is the cave depth from the surface. We
monstrous form may assume the solution in a then assume that the cave in S” has exactly the
non-stationary case, if we want to consider, for temperature of the rock at depth H in the system
instance, the cooling of a mountain during the S’ plus T0.
karst creation inside it? Let us consider now the three temperature
Nevertheless there is a simple way to lead fields. They are the solution of general Fourier
this problem to typical situations of engineering equation (Isachenko, 1969)
thermal exchange.
1 ∂T
Let us show how, considering three different ∇ ⋅ ∇T ( x, y , z ) =
systems, S, S’ and S”: a ∂t
1) The main system just now described, the That in our case, not time dependent, it
real case, which we call from hereafter S (Fig. reduces to the Laplace equation
8);
2) The system composed only of semi-infinite ∇ 2T = 0
undisturbed rock, without caves and with
It means that the T fields behave like a huge
external temperature equal to 0, which we call
class of phenomena for which the sum of the
from hereafter S’ (Fig. 9);
three spatial second derivatives is zero. The
3) The more complex system S”, (Fig. 10), functions that satisfy these conditions are called
composed by a cave at a particular temperature “Harmonic Functions”, and are among the most
T” buried in a semi-infinite rock, that releases important and studied functions in Physics
energy to the surface at temperature T0. Also the (Carslaw and Jaeger, 1959; Bejan, 1993;
rock at the infinite is assumed to have Balcerzak and Raynor, 1961; Nashchokin,
temperature T0. In this last ideal system S” there 1979). By the way, it would be possible to solve
is no geothermal flux at all. our fields using solutions given for different and
We are then ready to do the final step, well-studied problems like, for instance, the
assuming that T” in S” has a value given by: electric field due to particular charge
distributions, but we can do better for our
Fgt purpose.
T " = T0 + H [3.2]
Kr

Fig. 9. Qualitative view of the stationary thermal Fig. 10. Qualitative view of the stationary thermal
field due to the interaction between the geothermal field between a “hot conduit” at temperature T” and
flux and a conduit with no drainage (system S’). The the surface at T0 (system S”). The isothermal
isothermal surfaces are unaffected, and the cave surfaces are finite, closed and contain the conduit.
temperature is higher than T0. Heat diffusion problems like this one are quite usual
in thermal exchanges engineering.
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.12

The field structure around S, that we have K R ∇T = K R ∇T '−K R ∇T "

qualitatively shown (Fig. 8), is our unknown
term. These terms are now the energy fluxes that
flow through the systems S, S’ and S”, and then
The field structure of the second S’ it is r r r
obvious, it is composed by many horizontal F ( x, y , z ) = F ' ( x, y , z ) − F " ( x, y , z )
parallel lines (or planes, in 3-d) for T-field and
vertical vectors for thermal fluxes. But the energy flux in the system S’ is simply
The third system is the most exciting. There -Fgtk where k is the unit vector in the z-
are almost no lines all around (neither for T nor direction, and then
for thermal fluxes) unless in the region between r r r
the cave and the surface, because the F ( x, y , z ) = − Fgt k − F " ( x, y , z )
temperature difference between the “hot cave”
and the surface and also the “infinite”, drives an This means that we are able to calculate the
energy flux. From the other side, this situation is energy fluxes in the system S with vectorial
a very usual situation for engineering, because subtractions between the S” system, complex
the “hot cave” can be a tube transporting hot but already studied, and the trivial S’. If we
fluid buried in some engine... multiply this equation by the surface element dS
and integrate on a wide surface A that contains
Now we can do the last step: We state that the all the surroundings of the cave we have
unknown T field of the system S is given by
FA = − Fgt A − F " A
T ( x, y , z ) = T ' ( x, y , z ) − T " ( x, y , z ) [3.3]
The term F×A describes the flux outgoing
That is, we can subtract the (very complex, from the surface in presence of the cave, Fgt×A
but very common) T” field from the trivial T’ to the total flux if it would not be the cave, than
obtain our solution T. the energy flux captured by the cave is the
It is possible to prove this theorem in three difference between the two
steps:
Wcapt = Fcapt A = Fgt A − FA = [3.4]
1) The Laplace equation is linear, then if T1
and T2 are solutions, also (T1-T2) is a solution: = Fgt A − Fgt A + F " A = F " A
here in particular T’ and T” are solutions, then In this way the problem of energy
also T is; interception of a cold cave buried in an energy
2) The boundary conditions written above are flux is reduced to the energy transfer between a
satisfied by a T field given by (Eq 3.3); hot cave and the surface.
3) Then T is a solution of our equation with Now it is possible to study the T” solutions,
these particular boundary conditions, but the going to the heat transfer engineering to use its
solution is unique for the Uniqueness Theorem results.
for Harmonic Functions then, T is the solution...
The base of this proof is the linearity of the The shape factor
grad operator, which allows the first property. The thermal transfer engineering uses a very
But also the temporal derivative is linear, and effective approach to the problem of complex-
then we have another independent fundamental shape systems.
result: The T field may be calculated by this
way also for transient conditions, if we use the Let us return back to the fundamental
equivalent transient solutions for the T”. We are equation of conductivity, now written in three
not ready to use here this corollary, probably dimensions. The thermal flux through a surface
very important. The geothermal energy focusing element dA is given by
on caves dW = K R ∇T dA
It is interesting to look for other
consequences. Let us return to our [Eq. 3.3] to Where KR is the body (rock) conductivity.
apply the grad operator (that in equation will be We may consider two sources at definite
noted with ∇) and multiplying for the rock temperature T1 and T2. The thermal energy is
conductivity KR drained between the two by an intermediate
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.13

temperature field configuration, which depends T − T 

on the system shape in an extremely complex W = −K R A 2 1 
 ∆z 
manner.
Let us call Aiso the isothermal surfaces that We see that the integral behaves like the ratio
we can draw between the two sources; these between A and the sources distance, and then
surfaces can be infinite, also if neither of the the system geometry is included in this term
two bodies is infinite.
A ∂T *
The grad (T) must be perpendicular to these  =
 ∆z 
∫ Aiso ∂n
dA
surfaces because the thermal flux vector cannot
have any component along an isothermal The term in brackets is the ratio between the
surface, and then the thermal energy flows surface crossed by the thermal energy and the
normally to these surfaces. distance between the two sources. It is a
Let us call n the coordinate along the thermal “length” that characterises each system shape
flow path. The previous equation gives the that exchanges energy among two sources.
thermal flux as This scale-length is called “shape factor” in
∂T literature (Carslaw, 1959; Hahne, 1975;
dW = −K R dA Holman, 1996; Ozisik and Necati, 1993; Kays,
∂n 1966). We adopt unwillingly the usual notation,
The total flux between the two sources is then that uses the “S” for a length, but we shall write
given by the integral of thermal flux on any of SF, hoping to reduce (perhaps…) confusion with
these surfaces. It is not important to choose one the subscript “F”. Then
or the other, because we have assumed that the ∂T * [3.5]
only two “heat producer or destroyer” are these SF = ∫ Aiso ∂n
dA
two sources, and the energy must be conserved.
The thermal transfer through one of these The geothermal power absorbed by the cave
surfaces is therefore can be then written as
∂T
W = −K R ∫ Aiso ∂n
dA W = K R S F ∆T [3.6]

Now we have to study the way to use this

We define now a new dimensionless result.
temperature T* (that is in fact a relative
temperature variation in the path between the The “shape factor” calculation
two sources) as
We do not have to study in details the way to
T − T1 calculate the shape factors. Still in simple
T* =
T2 − T1 configuration the isothermal surface calculation
is very complex and the integration is in general
Where T1 and T2 are the sources extremely difficult, but the heat-exchanges
temperatures. Then it is possible to write literature contains many shape factors worked
out for the most common geometrical
∂T *
W = −K R (T2 − T1 ) ∫ dA configurations. Most of these results have been
Aiso ∂n obtained based on advanced analytical methods
(conformal mapping, superposition, special
In this way the sources temperatures are transforms, analogies with the electrical
analytically separated from the system potential studies and so on); it is not useful for
geometry, which now is completely included in us here to study these approaches.
the last integral, which is nevertheless extremely
Table 2 shows these shape factors in
complex also for trivial configuration.
interesting situation (Holman, 1996)
This equation has to be compared with the
We are going to use these shape factors, but
usual equation which describes the thermal
first it is necessary to answer an important
energy transfer between two sources separated
question: Are the shape factors linear? Let us
by a uniform distance ∆z through an area A discuss the question with an example.
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.14

TABLE 2
Semi-infinite medium with 4πR T1
isothermal surface and isothermal SF =
spherical cavity of radius R at depth 1 − R / 2H
H
H
T2
R
Semi-infinite medium with 4πR T1
isothermal surface and isothermal SF =
disc of radius R parallel to the π / 2 R − tan −1 (R / 2 H ) T2
surface at depth H H
R
Semi-infinite medium with 2πL T1
isothermal surface and isothermal SF =
cylindrical cavity of length L of cosh −1 (H / R ) T2
radius R, parallel to the surface at
H
R
depth H
L
Semi-infinite medium with 2π H T1
isothermal surface and isothermal SF =
cylindrical hole of radius R drilled to ln (2 H / R ) T2
a depth H normal to the surface.
H
R
Semi-infinite medium with 2π W T1
isothermal surface and isothermal SF =
plate (width W, long L, H>>W) ln (4 L / W ) T2
parallel to the surface at depth H H
W L

If we know the shape factors of two large that each temperature field is not affected
independent systems S1 and S2, say, one by the other. Only in this latter case the shape
composed by two cylinders and the other of a factor of the two conduits is the sum of the
sphere and a cylinder, can we consider a third shape factor of two single conduits, but in
system S3 (in this case, two cylinders and a general it is not so.
sphere) as composed by some “sum” of the two
firsts, and consider that its shape factor is given
4. The interception of geothermal flux by
by some “sum” of the two?
caves
Unfortunately the answer is: No, we cannot.
The temperature fields are linear and then the S3 The geothermal cross sections of caves
temperature field can really be calculated from
the S1 and S2, but it changes completely the It is now possible to apply the previous
equipotential surfaces on which the integration results to the problem of interaction between
is performed to “average” the flux in the caves and geothermal flux. It is stated above
integral [Eq. 3.5]. It is therefore necessary to that, to satisfy the boundary conditions, the T”
recalculate these surfaces and to repeat the temperature has to be
integration that will give a result that has no
Fgt
direct connection with the integrations of S1 and T " = T0 + H
S2 fields. Kr
This means that, for instance, the knowledge
As a consequence, the equation that gives the
of the shape factor of a conduit buried in a semi-
power intercepted
infinite medium tells us almost nothing about
the shape factor of two parallel conduits in the W = K r S F ∆T
same medium, unless their distance would be so
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.15

2π L 1900
SF = = = 250 m
−1  H  7.6
cosh  
 r 
With Fgt=2×106 J m2a-1 we have

W = Fgt SF H = 2 ×10 6 (250 × 500) = 2.5 ×1011 J a −1

Which is really a big power. It is possible to

study how much this figure changes with
conduit radius. Table 3 shows the energy
collected, by conduits of different sizes.
The result in the third column at first appears
surprising, because it shows that the variation of
the conduit size does not affect so much the
Fig. 11. A deep draining conduit can create a very intercepted power, but it is reasonable because
large geothermal shadow on the surface. Its scale the power is not absorbed by the conduit
size is not the conduit size, but the conduit depth. surface, but by the focusing effect of conduit on
the temperature field.
It is reduced to the very simple result The fourth column shows the surprisingly
high average thermal flux (note that are Watt
Fgt per square metre!) that enters through the
W = K R SF H = Fgt SF H
Kr conduit transverse surface. For comparison, the
[4.1]
Sun deposits on average 1.4 kWm-2 on the Earth
That is surface: The geothermal energy deposition on
small conduits is then of the same order! This
Intercepted flux = appears to be absolutely unbelievable, but is it
= (geothermal flux ) × (caveshape factor ) × (cavedepth ) true?
Table 3 shows that the scale-size of shape Roughly, the answer is that: Yes, it is true.
factor is some 1-10 times larger that the scale- But there are other important details to be taken
size of the underground structure. The last into account.
equation has to be compared with [Eq. 3.1].
The heating of water in deep conduits
Therefore, the effective area (we call it
“thermal cross section”) for thermal flux Let us discuss the effective water heating in
absorption of an underground conduit (Fig. 11) the focusing conduit, calling T its temperature at
is not its geometrical area but instead SFH, the springs.
roughly 10 times the cave size multiplied by its
depth, then it is always enormously greater than In the previous discussion, it was made the
the cave’s actual area! This amplification is due fundamental assumption that the water
to the “converging lens” effect created by the temperature T0 in the conduit does not change
presence of cold fluids in the cave that affects and really it is its low temperature that changes
the whole structure of the rock temperature the whole temperature field of surrounding rock.
field. This is equivalent to assume that the water flux
(or air flux in case of dry caves) is so large that
For example, let us estimate the geothermal the enthalpy intercepted by the conduit flows
power intercepted by a conduit at a depth away in the form of a small temperature
H=500 m with radius r=0.5 m and long L=300 increase of a very large amount of fluid, and
m. We may use the shape factor given by does not really affect the conduit temperature.
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.16

TABLE 3
Conduit radius r Shape factor SF Intercepted power Average Flux on
[m] [m] [J a-1] conduit [W m-2]
0.01 165 1.7×1011 900
0.1 210 2.1×1011 110
0.5 250 2.5×1011 25
2 310 3.1×1011 8
5 360 3.6×1011 3.8

We have then that the maximum energy water is minimum and the water flows out quite
extraction efficiency is obtained if the exit hot, at T”.
temperature T is almost equal to T0. Let us define the “critical fluid flux” Qc that
We have the opposite limit, if the warming is divides these two scenarios in a usually
so large that the water temperature T becomes idealised way. We look for a water flux Qc that
T”, the temperature of undisturbed rock. In this enters at temperature T0 and flows out at T” in
case the conduit becomes completely stationary conditions. The enthalpy subtracted
“transparent” to the geothermal flux, the to the system is
temperature field assumes a regular geothermal
gradient not affected by the cave presence and,
dE = C w Qc dt (T "−T0 )
as a consequence, just a little energy is If we admit that the system parameters do not
intercepted. The water into the conduit is then a depend on time, this enthalpy deficit must be
hot “mine water” and the classical, “wrong” given by the incoming geothermal flux Wdt.
solution [Eq. 3.1] becomes correct. Then using [Eq. 4.1] we have
The real cases are intermediate between these
two extremes, because in first approximation the Fgt SF H dt = C w Qc dt (T "−T0 )
water really warms but, as a consequence, its
capability to intercept geothermal energy is But T” is given by
reduced, because the temperature difference
Fgt
between the rock and the water becomes T " = T0 + H
smaller. Really the scenario is still more KR
complex, because the water temperature
increases along its path underground: It enters And we obtain
cool, very efficient in geothermal energy
[kg s−1 ]
dM K R SF [4.2]
focussing, but as long as it warms downstream Qc = =
dt Cw
its capability to intercept the geothermal flux
decrease. The non-linearity of SF forbids correct Then the critical flux is simply the conduit
analytical solutions, but we can make some shape factor “scaled” by a dimensional term
other step. (rock conductivity divided by the thermal
capacity of flowing fluid). In terms of volume
The critical shielding discharge flux
It is possible to estimate the water heating
[m 3s − 1 ]
dV Qc K R SF
along a deep conduit and its final temperature as = =
a function of conduit parameters. dt ρ w C w ρ w
We have just seen that two extreme scenarios From another point of view, Qc can be
are possible. If the fluid flux is very large the considered the critical flux below which it is
rock temperature field is completely changed, possible to consider that the rock temperature
the geothermal energy flow interception is field is undisturbed. Or, from still another point
maximum and the water flows out at T0. At the of view, we can be sure that a water flux much
other extreme, if the water flow is very small, larger than Qc perturbs the rock.
the rock temperature field is completely
undisturbed, the geothermal energy flow to the
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.17

The solution [Eq. 4.2], which is surprisingly We are going to improve this point of view
simple in comparison with the bloody analytical and discover that Qc has another, still more
difficulties of the problem, is able to give also important, meaning.
the velocity that water must have to satisfy the
equation. If the conduit radius is r, then Geothermal power intercepted
1 K S 1 We can calculate the effective cave
v = Qc = R F temperature T at the equilibrium and solve the
π r ρw
2
Cw π r 2 ρ w
inverse problem, the estimation of flowing
In the case of a water draining conduit in depth of hot spring waters.
limestone we have The final system temperature must lie
between T0 (near it for high water discharges)
2.3
Qc,water = SF = 5.5 ×10 −4 SF [kg s−1 ] and T” (near it for low discharges). Let us call
4.2 ×10 3 [4.3] WM the maximum power that it can be
subtracted by our system
If the conduit drains air, the critical flux
becomes W M = Fgt SF H = 0.06 SF H [W]
[4.4]
2.3
S = 2.3×10 −4 SF [kg s−1 ]
Which is valuable for very large water flows
Qc,air =
10 3 F and outflowing temperature around T0.
If we call T the real (unknown) conduit
The shape factor is in general comparable
temperature at the equilibrium, the residual
with the conduit length, usually some 103-104
outgoing upward flux is not zero, because it
m, therefore the water flux able to create
“filters” an energy given by
changes in the rock temperature fields is in
general quite small. Neither the air flux Wup = K R S F (T − T0 )
requested to cool the rock it is too large, in
absolute, but its small thermal capacity and its The difference between WM and Wup is the
very small density cause a large volume flow net flux entering inside the cave from below. If
request. the conduit is at temperature T, the energy
It is nevertheless necessary to remember that conservation states
temperature field changes are only possible if
the fluid fluxes have had sufficient time (more Win = WM − Wup = Fgt S F H − K R S F (T − T0 )
than ∆teq defined in [Eq. 1.1]) to converge at the
stationary (equilibrium) state. But from [Eq. 3.2] we have
The critical flux Qc for air or water are Fgt H
extremely important for another reason: They KR =
are the air or water fluxes able to shadow the T "−T0
upper rock from the geothermal flux, forcing its
Then
temperature near to the average yearly
temperature of external atmosphere, Tave.  T − T0   T "−T 
Finally, it is important to note that this work Win = Fgt SF H 1−  = WM  
 T "−T0   T "−T0  [4.5]
was originally a chapter of an underground
climate physics book. Its purpose was to
If T=T0 the equation reduces to Win=WM
estimate when we could consider a cave as
(which describes the “system S” situation), and
“shielded” from the upward geothermal energy
if T=T” the term Win vanishes, as expected.
flux. If the flux is much larger than Qc we are
sure that the mountain parts (and caves) above Therefore, the geothermal flux intercepted by
the conduit are shielded; if the flux is well the cave is reduced as long as its final
below the critical value the conduit is temperature T increases: The cave is becoming
“transparent” to the geothermal energy and we “transparent”, and this equation describes its
have to include also its contribution to analyse “fading” inside the temperature field.
the underground climate in the rock above the
conduit.
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.18

Temperature and deposited power versus the ratio between this difference and the
discharge theoretical, maximum difference T”-T0. Then
Usually we have very small possibilities to T "+qT0  1 
establish a natural conduit discharge. We deal T − T0 = − T0 = (T "−T0 ) 
with actual discharge Q, and we can only try to 1+ q  1+ q 
understand how this water (or air) flow is
And calling the “excess temperature ratio” of
affected by the geothermal flux.
groundwater, that is the amount of actual
Let us then return to the application of the heating in comparison with the maximum
First Principle [Eq. 4.5]; if we call Q the total attainable, we have
actual water flux in the system, it gives T − T0 1
rT = = [4.10]
 T "−T  T"−T0 1+ q
WM   = C w Q (T − T0 ) [4.6]
 T "−T0  With the assumption [Eq. 4.9] we can rewrite
Where we have assumed that the whole the [Eq. 4.5]
entering energy flux goes to warm the water
flux. This statement is true if we admit that the  T "−T  q
Win = WM   = WM [4.11]
thermal energy is released on the whole system,  T "−T0  1 + q
i.e. it is true if the temperature T, which is
actually the output water temperature, can be We have then two fundamental equations, [Eq.
used to describe the whole system, also in its 4.10] and [Eq. 4.11], which connect the internal
further parts. It is a strong, and in general false, drainage Q to the outgoing temperature and to
hypothesis, but it is better to assume it as true the intercepted geothermal flux. The two
and only afterwards have a look on what graphics (Fig. 12) describe the behaviour of the
happens in more real situations. out-flowing water temperature T and of water
With the previous results and trivial absorbed energy as a function of discharge,
calculations, we have obviously measured in function of our nice
scale-discharge Qc (it can be adapted to air flow
WM  T "−T  with trivial changes).
Q=   =
C w (T − T0 )  T "−T0  We have previously discussed the Qc as the
[4.7] “shielding flux” and WM as the “maximal
WM  T "−T   T "−T  intercepted flux”. Now we see that they are
=   = Qc  
C w (T "−T0 )  T − T0   T − T0  mainly the natural scales of fluid flow and of
geothermal power flux intercepted, exactly as
happens with the Similarity Numbers, always a
This important equation relates the actual flux ratio between a parameter and a scale-value for
Q to the output temperature T in terms of the it.
critical shielding flux Qc and the surrounding
temperatures. It is very easy to solve it to obtain
Q
T "+T0
Qc T "+ qT0
T= = [4.8]
1+
Q 1+ q
Qc

Where we have called

Q
q= [4.9]
Qc
Fig. 12. Water temperature increase and intercepted
The excess temperature above the “external power by a deep conduit versus water discharge.
average” T0 is then T-T0, but its natural scale is
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.19

Really, the apparition of a natural scale-value  Q 

of the main variables, the fluid flow and the H = H 0 1 + −4

geothermal power, may suggest to use this  5.5 × 10 S F 
typical thermo- or fluid-physicists approach, This solution shows that the evaluations made
calling the Q/Qc and W/WM terms as Someone with the Desio formula [Eq. 4.12] are deeply
Number, to create two underground companions underestimated, unless for very small
to Nusselt, Reynolds, Peclet, Froude and so on discharges. In fact we can write
Numbers.
Nevertheless we do not like this way to H = H 0 (1+ q) [4.13]
describe physical processes, preferring It is then possible to see that the q-number is
(aesthetically, because mathematically it is essentially the “amplification” term of estimated
absolutely equivalent…) to work with the scale depth H0.
dimensions. Therefore we do not propose this
Similarity approach although, if a name must be The main difficulty in these formulas it is the
given to the geothermal “heat” term W/WM, we estimation of critical discharge Qc, which
strongly suggest “Alighieri Number”... requires the knowledge of the conduit shape
factor, in generally unknown. In case studies it
The inverse problem is necessary to take into account the geological
context to estimate the probable conduit shape
It is obviously of main interest the problem of in order to calculate the critical discharge Qc.
deducing the provenance depth of a hot spring.
Let us do an example. In many case, for
The traditional Desio formula (Celico, 1986)
instance, we can assume a “U” shape for the
assumes essentially that a water flux at
whole drainage system. A similar conduit can
temperature Tw comes from a depth H0 at which
simply be approximated with a conduit of length
the temperature of undisturbed rock is Tw, that
L at depth H, because its two vertical branches
is
do not matter for the shape factor, being merged
in a rock shielded by the deep drainage. Then
H0 =
KR
(Tw − T0 ) [4.12]
Fgt 2πL
SF = [4.14]
This is true for “mine waters”, deposits of cosh −1 (H / R )
resident waters, which have essentially the T”
The Table 4 gives the value of inverse
temperature, but we have seen that this
hyperbolic cosine term for typical conduit
assumption is in general false, because a natural
radius and depths.
water flux had surely succeeded to disturb the
whole rock temperature field. We have then to
TABLE 4
use [Eq. 4.8] and [Eq. 4.2] to write
H↓ R→ 0.1 m 1m 10 m 100 m
−1
  200 8.3 6.0 3.7 1.3
Tw − T0 = (T "−T0 )1 +
Q
 =
500 9.2 6.9 4.6 2.3
 Qc 
−1 2000 10.6 8.3 6.0 3.7
Fgt  Q 
= H 1 + 
It is possible to see that, unless the nearby
K R  5.5 × 10 − 4 S F  “pathological” cases of the last column, the
That correlates the temperature increment denominator in [Eq. 4.14] is not far from 2π;
above the local temperature to the unknown therefore for this conduit it is possible to assume
depth H, to the discharge Q and the conduit a shape factor roughly equal to its length
shape factor SF. Then
SF ≈ L
K  Q 
H = R (Tw − T0 )1+ −4  The formula that estimates the water
Fgt  5.5 ×10 SF  provenance depth can then be approximated as

Comparing with [Eq. 4.12] we obtain the  Q 

H ≈ H 0 1 +  [4.15]
solution for thermally disturbed rock (Q average −4
 5.5 × 10 L 
yearly discharge)
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.20

In general a typical deep circuit has a very It is possible to perform a last step,
large size, many kilometres. If we call L* its considering a long conduit L, along which the
length expressed in kilometres we have, for the water is heated. The term T is now the
circuit described above temperature at the length x, in a section dx with
shape factor sF. The thermal flow and the
Qc ≈ 0.5 L * temperature increase in that section is then
given by [Eq. 4.5]
And finally
dW = dWM − dWup =
 2Q 
H ≈ H 0 1 +  [4.16] Fgt H
 L*  = Fgt s F H − s F (T − T0 ) = C w Q dT
T "−T0
It is a simple formula that estimates the
drainage depth as a function of discharge. And then
As an example, let us consider a spring with Fgt S F H  T "−T  sF
an average discharge Q=20 kg s-1, and a   = dT
temperature 20°C above the local average. The CwQ  T "−T0  S F
depth estimation [4.12] gives H0=830 m.
This is an equation that could solve the
If the hydro-geological context suggests a problem, if we would able to integrate the left
circuit length L*=20 km, the critical discharge part, but unfortunately this is not possible.
is Qc=11 kg s-1, and then the q term is almost
equal to 2. Our formula estimates then a depth We have written sF and not dSF (as would be
H of 2.5 km for the circuit, much more than natural) because it is not possible to pass from
supposed… the equation that gives SF as a function of L, to
the contribution of a part dx of L to SF. At the
With this value we can return to the end of the third chapter, we have noted that SF
geological context and, if we have other does not linearly depend on each dx part,
information (like the circulation time), to more because it comes from an average on the whole
accurate estimations of the ratio between the space and system, and it is not possible to
depth and the radius of a deep drainage system. consider it as the result of an integration on
A better estimation of provenance depth with some dx.
[Eq. 4.13] and [Eq. 4.14] is then possible.
For instance, the contribution of the dx at one
It is nevertheless better to remember that we conduit edge covers a cone above it, and the
are working with the assumption that the system local sF is like that given by a small sphere,
has attained stationary condition; the above whereas the dx in the middle of conduit gives a
formula is then correct for water fluxes that very small contribution.
persist from very long times, much longer than
We can nevertheless integrate the last
∆teq of [Eq. 1.1].
equation assuming the false approximations
sF=dSF. It makes no analytical sense but it
Temperature changes into the system probably introduces an error smaller than the
It is possible to perform a last “calibration”. assumption of uniform system temperature.
We have already noted that we are assuming Then for a conduit buried in a semi-infinite
that it is possible to define univocally a system medium at depth H
temperature T, but this is not always true. It is
Fgt S F H  T "−T  dx
possible to speak univocally of “system   = dT
temperature” if each part of the system is CwQ  T "−T0  L
uniformly heated by the energy flux, for
instance if the head and the tail of water flux are Using Qc becomes
mixed (for instance, when the water enters in a
spherical deposit). Nevertheless this is not the 1 Fgt H  T"−T  dx Qc dx
  = (T "−T ) = dT
usual situation because in a real conduit the q K r  T "−T0  L Q L
water enters with a temperature T0 and flows
warming up to the final temperature T. And integrating on x from 0 to L and on T
from T0 to T we have
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.21

 Q  water flows that have attained a steady state

T"−T = (T "−T0 )exp− c  situation with rock.
 Q
Returning to the previous example, of a
Adding and subtracting T0 from the equation spring with Q=20 kgs-1, and a temperature 20 °C
left side, we have above the local average, H0=830 m, Qc=11 kg s-
1
, and q=2.
 1
T − T0 = (T "−T0 ) − (T "−T0 )exp− The corrective term to be applied to H0 with
 q [4.17] [Eq. 4.16] it is a factor 3, but now [Eq. 4.19]
gives a factor 2.54.
This has to be compared with [Eq. 4.10]. It is
Really the temperature variability along the
a simple and nice equation that describes the
conduit gives [Eq. 4.19] a final heating at a
water heating during a flow.
depth smaller than in the case of a “global”
The Qc term has returned, and continues to be heating [Eq. 4.13], but the difference does not
the scale discharge of deep conduits. If the appear as too significant if compared with the
effective discharge Q is large (in comparison intrinsic uncertainties of such problems.
with Qc) the water temperature at the outflow is
near the T0, and if q is near 0 the T=T”. Steady State Geothermal Power Plant
The shape factor has disappeared, because
In the previous chapters, it has been shown
this equation is valuable everywhere the (strong
implicitly a way to extract power from
and false) sF=dSF approximation is valuable. underground, using a deep conduit that focuses
Nevertheless we can suppose that it is on itself large amounts of geothermal energy.
reasonable model, and we guess that the last This is deeply different from the usual
equation gives a fair approximation of natural Geothermal Power Plants, which extracts
heating processes along a conduit. energy (or, better, are believed to extract,
It is useful to invert again the problem to because cool water always focus on itself the
obtain the estimated depth crossed by water of a temperature field) from hot rock, directly
spring at temperature Tw. We have cooling it. In principle, when the rock is cooled
the power plant stops its work.
Fgt   1  Here we have shown that the deep cooling
Tw − T0 = H 1 − exp − 
Kr   q  effect acts as an energy attractor on the cooled
rock, and then that a power plant working in
And with simple passages such way, it will never end its fuel.
We want here to make the next step, looking
−1
K (T − T0 )   1  for its “constructive” efficiency.
H= r w 1 − exp −  =
Consider a fluid that transfers energy QH
Fgt   q 
from a “hot” source at TH to a colder source at
H0 [4.18] TL. Does this process produce work? If the
=
 1 energy transfer is made with “special” systems
1 − exp −  it does: They are called “thermal engines” and
 q use the energy flow from TH to TL to produce
That corresponds to [Eq. 4.13], and reduces work. A power plant is said to “produce”
to it for q~0 and q>>1. With the same energy, but this is trivially false because the
assumption of [Eq. 4.16] on shape factor, and energy cannot be neither produced nor
L* in [km], it gives destroyed: It stores energy at a very low entropy
(“work”) in an entropy flow from a high
H0 temperature (low entropy) to a low (high
H=
 L*  entropy).
1 − exp −  [4.19] The Second Principle of Thermodynamics
 2Q 
states that the efficiency -that is the ratio
between the work rate given L and the heat rate
That can be considered a reasonable formula
absorbed W- of a reversible thermal engine
to estimate the water provenance depth, for
working between the two sources is
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.22

LMax TH − TL T It is easy to calculate the value qM for which

= = 1− L the entropy flow attains its maximum
QH TH TH
T "−T0 T"
But we have that q M = γ +1 = +1 =
T0 T0
QH  TL 
∆S = 1 −  That we can substitute in the previous
TL  TH  equation to obtain the maximum of entropy
flow. In natural cases the term T0 is some 280
And then
K, the T” some 350K, then the ratio is slightly
LMax = TL ∆S more than 1, and then

That is, the maximum work available is given T"−T0 1 F 

by the product of the temperature of cold source qM = 1+ ≈ 1+  gt H
T0 2  K rT0 
and entropy change during cooling, which is the
Free Energy variation in the transformation. Which gives, erasing the second-order terms,
In the case of interaction cave-geothermal the maximum power production of this
field, the cave acts as a thermal sink, in two geothermal power plant
phases. In the first step it intercepts a flux of
geothermal energy W (low entropy, temperature  ∆L  1  T"−T0 
T”) from downward, which results in a water   = T0 ∆SMax = W M  
temperature increase from T0 to T. In the second  ∆t  Max 4  T0 
step the energy it is released as “disordered
energy” to the atmosphere (high entropy, T0) at It is necessary to emphasize the difference
the spring. For instance, considering the deep between the subtracted power WM and the
conduit as a geothermal power plant, we have maximum available work (or power) LMax. The
first is interesting to make some use that
that its entropy production per time unit ∆t -we
requires enthalpy at constant temperature, as it
are dealing with discharge Q- is
is the case of ice melting or water evaporation,
 W W  q  1 1  for which WM, not LMax, is used. But to create
∆S = − + ∆t = SF Fgt H   − + ∆t structures we need “work” also in the physical
 T T0  1+ q   T T0  sense: Order, available work.
Therefore, the LMax terms in each water
Calling heating and rock cooling processes are directly
connected with the entropy rate at disposal for
T "−T0 Fgt H constructive processes, that is, they may appear
γ = =
T0 K rT0 as the building rate of ordered structures, like
conduit networks.
With some work and using [Eq. 6.1]

∆S W M  q  1  Geothermics and phreatic conduit genesis

= γ  
∆t T0 1+ q  γ +1+ q  We have observed above that the initial
purpose of this work was to show that the
This shows that the entropy production goes geothermal energy flux could not participate in
to zero for q=0 and q=∞, because if the water the characterisation of cave climate and then,
flow is very small the fluid final temperature is for instance, to speleogenesis (Badino, 2005).
quite high but the total energy removed is very As the reader has seen, we are showing
small; from the other side, if a lot of water flows exactly the contrary.
into the conduit, its final temperature is
essentially T0, then the entropy is able to flow Here we are going to make the last step
between the rock and the water, but it is not giving some ideas about the geothermal role in
finally transferred to the atmosphere and to an the genesis of phreatic conduits and in general
external “final user”. of underwater drainage networks.
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.23

During deep flow the underground waters are carbonate dissolution releases ions that are in
warmed from their starting temperature T0 up to part the same already present in water.
the final T, that has a theoretical maximum at This complex system forces to find the
T”, as shown above. solution of many different equations describing
What is the typical temperature increase? We equilibrium kinetics, charge and mass
have seen that a conduit is able to shield the conservation. General solution charts are given
geothermal flux like a plane watertable does, in (Ford and Williams, 1989); they show the
then the scale temperature increase is around the saturation values at the equilibrium for various
value given by [Eq. 2.1] initial partial carbon dioxide pressures. It results
that in open systems (with release of carbon
0.06 500
∆Tgt = = [°C] dioxide excesses) the warming of a calcite
4.2 ×10 (P − Pout ) P *
3
saturated water gives, without exceptions, a
super-saturation and then provokes a calcite, or
Therefore, in real cases it is in the range aragonite, deposition.
between 0.2 and 3 °C, a temperature drop that
In a closed system this behaviour changes in
water gains during flow between the cave
a complex way. Generally a super-saturation is
bottom and the springs. Now we know that this
produced, but if the initial CO2 partial pressure
temperature change happens also along
is below 200 Pa (0.002 atm) and the temperature
conduits, not only in the “watertable”, and that
is below 30 °C, a calcite under-saturation
the power release is concentrated on the conduit
appears as result of water heating, as larger as
surface walls. How does this warming affect the
colder is the water.
water chemistry?
The typical carbon dioxide partial pressure in
We can outline the chemical behaviour of
free atmosphere is 3.5×10-4 atm, then at 10 °C
water saturated of calcium carbonate entering in
the calcite equivalent content at the saturation is
flooded conduits.
around 12 mg l-1. A water temperature increase
It is well known that the carbonate of 1 °C result in a saturation value of 0.02-0.04
dissolution in water is very complex (Snoeyink, mg larger, i.e. with a flux of 1 m3s-1 it gives
1980), (Ford and Williams, 1989), (Dreybrodt, around 103 kg of dissolved rock per year.
2000) because the equilibrium state results from
It is a small figure but it suggests that further
the combined equilibriums of different,
studies are necessary to a more complete
interconnected reactions, which depends on
understanding of saturation conditions as a
temperature, local pH and the presence of other
function of temperature, of chemically complex
dissolved salts with common ions.
waters in a closed system.
In the simplest case, the first equilibrium
In any case indirect evidences of
reaction gives the amount of dissolved carbon
effectiveness of speleogenetic processes
dioxide, for which in usual conditions the Henry
induced by geothermal heating in phreatic
Law holds, stating that the dissolved gas
conduits can be found, because if these
decreases with temperature and it is
processes are possible, they have to affect the
proportional to its partial pressure above the
network morphologies:
water surface. Therefore, its quantity depends
also on the kinetics of gas transport until the 1) the geothermal energy is released only in
surface, if it does exist. the lowest conduit walls, then the dissolving
characteristics have to depend on the rock
The other reactions, which involve only water
surface orientation;
and carbon dioxide, are the dissociation of
carbonic acid in calcium bicarbonate and H+, the 2) a deep conduit shadows completely the
dissociation of bicarbonate and the equilibrium upper rock, then the formation of a conduit that
H+ and OH- in water. These dissociations tend cross the rock above another is hampered, and
to increase with temperature thanks to the this affects the whole drainage conduit structure.
increase of available energy. Similar processes can probably play a part
The last main reaction describes the also in the deep drainage network formation in
equilibrium between the calcium carbonate and glaciers (Badino, 2002), but either ice or
water enriched with carbon dioxide. The limestone, a lot of work has still to be done for a
better understanding of geothermal role in karst.
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.24

Conclusions Bejan A. 1993. Heat Transfer. John Wiley and

Sons.
The estimations of temperature fields inside
mountains are important for speleogenesis and Benderitter Y., Roy B. and Tabbagh A. 1993
for underground climate studies, but also for Flow Characterization through Heat
many cases which require an energy balance on Transfer: Evidence in a Carbonate Fractured
a sub-geological time-scale, like glacier Medium. Wat. Res. Res. 29, 11: 3741-3747.
stability, geothermal spring studies, deep hydro- Bohren C. and Albrecht B. 1998. Atmospheric
geological analysis, tunnel drilling and so on. Thermodynamics. Oxford Un. Press, 402 pp.
We have shown that these potentially Bogli A. 1980. Karst Hydrology and Physical
cumbersome modelling can be reduced to Speleology. Springer-Verlag.
simple calculations that allow quite accurate Carslaw H. and Jaeger J. 1959. Conduction of
estimations of energy absorbed by deep Heat in Solids, Oxford-Clarendon Press.
structures and of provenance depths of
Catalano P. 1993. Laboratori Sotterranei:
geothermal waters.
Relazione Geologica. Unpublished, INFN.
These results could also be applied for remote
Celico P. 1986. Prospezioni Idrogeologiche.
sensing of deep drainage structures and for
Liguori.
construction of inexhaustible geothermal power
plants, but at present these applications appear Dreybrodt W. 2000. Equilibrium Chemistry of
to meet insurmountable practical difficulties. Karst Water in Limestone Terranes. In
Speleogenesis: Evolution of karst aquifers,
Klimchouk, A., Ford, D.C., Palmer, A.N.,
References and Dreybrodt, W. (Eds.), Nat. Speleol. Soc.,
USA: 130-135.
AA. VV., 2004. Under the Desert: the
Mysterious Water of Cuatro Ciénegas. La Ford D. and Williams P. 1989. Karst
Venta. Geomorphology and Hydrology, Unwin
Hyman.
Badino G. 1995. Fisica del Clima Sotterraneo.
Memorie IIS 7, II. Goy L., Fabre D. and Menard G. 1996.
Modelling of Rock Temperatures for Deep
Badino G. 2000. I Gradienti di Temperatura nei
Alpine Tunnel Projects. Rock Mech. Rock
Monti, un Indicatore Esplorativo. Talp-FST
Engng 29.
21, 72-80
Guichonnet P. 1967. Il Traforo del Monte
Badino G. 2002. The Glacial Karst. Proc. of V
Bianco, Mondadori.
Int. Symp. on Glacier Caves and Cryokarst
2000, Nimbus 23, VII, 1/2002. Hahne E. and Grigull U. 1975. Formfaktor und
Formwiderstand der Stationaren
Badino G. 2004. Cave Temperatures and Global
Mehrdimensionalen Warmeleitung, Int. J.
Climatic Change. Int. J. Speleol. 33.
Heat Mass Transfer 18: 751-767
Badino G. and Forti, P. 2005. L’eccezionale
Holman J. 1996. Heat Transfer, MacGraw-Hill.
Ambiente della Cueva de los Cristales. Proc.
“Le Grotte di Miniera”, Memorie IIS XVII, Isachenko V., Osipova V. and Sukomel A.
II. 1969. Heat Transfer, MIR.
Badino G. 2005. Clouds in Caves. Jeannin P., Liedl R. and Sauter M. 1997. Some
Speleogenesis and Evolution of Karst Concepts about Heat Transfer in Karstic
Aquifers 2 (2), www.speleogenesis.info. Systems. 195-198.
Badino G. 2005. Nuovo Collegamento Kays W. 1966. Convective Heat and Mass
Ferroviario Torino-Lione, Temperature del Transfer. McGraw-Hill.
Tunnel di Bussoleno, Geodata, Torino, Lee et al. 1966. Heat Flow and Volcanic
unpubl. Temperatures. Handbook of Physical
Balcerzak M. and Raynor S. 1961. Steady State Constants, The Geological Society of
Temperature Distribution and Heat Flow in America.
Prismatic Bars”, Int. J. Heat Mass Transfer Koenigsberger J. and Thoma E. 1906. Uber die
3: 113-125, Beeinflussung der geothermischen
 Giovanni Badino / Speleogenesis and Evolution of Karst Aquifers, 2005, 3 (2), p.25

Tiefenstufe durch Berge und Taler. Eclog. Ozisik M. and Necati. 1993. Heat
Geol. Helv. IX (1). Conduction,.Wiley InterScience.
Laidler K. and Meiser J. 1995 Physical Schoeller H. 1962. Les eaux souterraines.
Chemistry. Houghton-Mifflin Co. Masson.
Lismonde B. 2002. Aérologie des Systèmes Snoeyink V.and Jenkins D. 1980. Water
Karstiques. CDS Isére. Chemistry. John Wiley & Sons.
Luetscher M. and Jeannin P.-Y. 2004. Szechy K. 1973. The Art of Tunnelling,
Temperature distribution in karst systems: Akademiai Kiado, Budapest.
the role of air and water fluxes. U.S. Bureau of Mines. 1996. Dictionary of
Speleogenesis and Evolution of Karst Mining, Mineral, and Related Terms, CD-
Aquifers 2 (2), (from Terra Nova 16):, 344– ROM.
350
Verhoogen J. 1956. Temperatures within the
Nashchokin V. 1979. Engineering Earth. In Physics and Chemistry of the
Thermodynamics and Heat Transfer. Mir, Earth, Pergamon Press, I.
573 pp.