Douce A.P. Thermodynamics of the Earth and Planets

Подождите немного. Документ загружается.

128 Energy transfer processes in planetary bodies

x+ x

x, y, z

+ y

, ,C

Fig. 3.3

Inﬁnitesimal volume element in which diffusive heat transport is occurring.

We use partial derivatives because temperature is in general a function of the three spatial

coordinates and of time. The temperature gradient in equation (3.8) is evaluated at x, which

is the coordinate of the element’s left face, but the temperature gradient is not necessarily

the same on the element’s right face. Given that the volume element has infinitesimal

dimensions, the heat flow across the right face can be written as:

−kδyδz

∂T

∂x



x+δx

=−kδyδz



∂T

∂x



∂

∂x



∂T

∂x



δx



=−kδyδz



∂T

∂x



∂

∂x

δx



. (3.9)

The net rate of change of enthalpy in the infinitesimal volume element due to heat flow in

the x direction is given by expression (3.8) minus equation (3.9) (i.e. “heat entering” minus

“heat leaving”):

kδxδyδz

∂

∂x

. (3.10)

Similar arguments for heat flow along the other two orthogonal directions show that the rate

at which the enthalpy of the infinitesimal volume element changes due only to the effect of

heat flow is:

kδxδyδz



∂

∂x

∂

∂y

∂

∂z



=kδxδyδz∇

T , (3.11)

where the operator ∇

, called the Laplacian, is shorthand for the sum of the second partial

derivatives of a function relative to each of the independent spatial variables (note that the

Laplacian is a scalar).

129 3.2 Heat transport by diffusion

If the material that we are considering contains heat sources or sinks, then their contribu-

tion must be added to equation (3.11) in order to obtain the total rate of change of enthalpy.

Let the amount of heat added per unit mass and per unit time be α. For a heat source it is

α>0, whereas for a heat sink it is α<0. At this point we do not care what the sources or

sinks are, or how many of them are present – the value of α is simply the net amount of

heat added per unit mass and per unit time, from all sources and sinks combined. The total

rate of change of enthalpy of the volume element, from the combined effects of heat flow

and heat sources and sinks, is, then:

=kδxδyδz∇

T +αρδxδyδz, (3.12)

where ρ is the density of the material (recall that α is rate of heat supply per unit mass).

The rate of change of enthalpy is also given by:

=ρδxδyδzc

∂T

∂t

. (3.13)

Equating (3.12) and (3.13) and simplifying we arrive at the heat diffusion equation:

∂T

∂t

ρc

∇

T +

(3.14)

or:

∂T

∂t

=κ∇

T +

, (3.15)

where κ = k/(ρc

) is another material property called the heat diffusivity, with dimension

[L]

[T]

−1

, e.g., units of m

−1

. Equation (3.15) is called the diffusion equation and its

mathematical form applies to transport of any quantity that is described by equation (3.4).

Momentum diffusivity is described by the kinematic viscosity ν given by the ratio ν =µ/ρ.

Chemical diffusivity will be discussed in Chapter 12.

The diffusion equation is an example of a partial differential equation. As all differential

equations, it does not tell us directly how to calculate a quantity. Rather, it gives us a

mathematical law that describes how the quantity of interest varies. In this case, all that

equation (3.15) says is that the rate of change of temperature is a linear function of the

curvature (second derivative) of the temperature field. Starting from this information we

can construct a precise mathematical description of the spatial and temporal distribution of

temperature. At this point you should study Box 3.1 in depth, as it is fundamental for much

that follows.

Box 3.1

Solution of the heat diffusion equation in a semi-inﬁnite half space

The heat diffusion equation is an example of a partial differential equation, as it relates the partial derivatives

of temperature relative totwo independent variables, space and time. Ingeneral, solving a partialdifferential

equation requires that we begin by specifying the geometry of the problem and a set of initial and boundary

conditions. The geometry of the problem is a precise description of the nature and size of the region of space

130 Energy transfer processes in planetary bodies

Box 3.1

Continued

over which we want to solve the equation, in this case, the region of space for which we wish to know the

temperature distribution at any instant of time. The initial condition is the value that the variable of interest

(in this case, temperature) takes at some instant chosen as time =0 and at every point in the spatial region

of interest. The boundary conditions are the values that the variable (or some of its derivatives) take at some

speciﬁed points in space (generally, along the boundaries of the region of interest) for all instants of time.

Here we will solve the heat diffusion equation in one dimension over a semi-inﬁnite half space, which is a

region of space that has only one boundary. This means that the temperature at a point inﬁnitely far away

from the boundary does not change with time. This particular solution is fairly straightforward and leads

naturally to the meaning of the characteristic length and time scales of heat diffusion. It is also an important

step in understanding heat transfer in terrestrial planets and in inferring their internal structures.

We seek a solution to the heat diffusion equation (3.15) in one dimension and without a heat source:

∂T

∂t

=κ

∂

∂x

. (3.1.1)

A solution to this equation is a function T = T(x,t). We imagine a semi-inﬁnite half space at a uniform

temperature and we impose an instantaneous temperature perturbation on the space’s only boundary. We

seekthefunctionthatdescribesthepropagationofthisperturbationinspaceandtime.Weneedtospecifythe

appropriate initial and boundary conditions. This is illustrated in Fig. 3.4. The initial condition (t

in Fig. 3.4)is

a semi-inﬁnite half space at a uniform temperature T

, except at the boundary (x=0) where the temperature

is T

. This is the temperature perturbation that we impose at time t

.IfT

< T

then the half-space cools by

losing heat across the boundary, as shown in Fig. 3.4. The temperature at the boundary is held constant at T

for all times – this is the boundary condition. Our initial and boundary conditions are, respectively:

(

x > 0, 0

)

(

0,t

)

. (3.1.2)

As the half space loses heat its temperature distribution changes as shown in Figure 3.4, where t

, t

are

three progressively longer intervals of time.

We wish to ﬁnd the function that describes temperature distribution in space and time, i.e. the function

that describes the curves in Figure 3.4. These curves can be thought of as sections of a surface extending

perpendicular to the page, where the third coordinate represents time. One technique for solving partial

differential equations (PDEs) is to convert them into ordinary differential equations (ODEs). These are

equations that involve derivatives (possibly of different degrees) relative to a single independent variable,

and are generally easier to solve. One way of accomplishing this is by combining the two independent

variables in the PDE into a single one. Because diffusivity has dimension of [L]

[T]

−1

we can deﬁne a

non-dimensional variable ζ (also called a similarity variable) by combining x and t as follows:

ζ =

√

κt

(3.1.3)

(the factor 2 in the denominator has no physical signiﬁcance, and is justiﬁed a posteriori because it simpliﬁes

the algebraic manipulations). Besides simplifying the solution of the partial differential equation, the

similarity variable gives us the scaling factor for length vs. time. Once we know what the temperature is at

any given ζ then we can immediately calculate all of the x–tcombinations at which temperature has

131 3.2 Heat transport by diffusion

Box 3.1

Continued

Fig. 3.4 Temporal evolution of the diffusive temperature distribution in a semi-inﬁnite half space whose boundary is kept at

a constant temperature T

< T

(the temperature at inﬁnity). At the initial time (t

=0) the half space has uniform

temperature T

, except at the boundary, where T =T

. The curves for longer intervals, t

> t

, are

differently-scaled versions of the error function (see text).

this same value. If we seek a function for temperature in terms of this non-dimensional variable, however,

then we should also express temperature as a non-dimensional variable. This is accomplished by deﬁning

the variable θ:

θ =

T −T

−T

, (3.1.4)

which scales temperature to the temperature range characteristic of the problem. In terms of the two

non-dimensional variables, our initial and boundary conditions become:

θ = 1 as ζ →∞

θ = 0 at ζ = 0. (3.1.5)

We now use the chain rule to ﬁnd expressions for ∂T/∂t and ∂

T/∂x

in terms of dθ/dζ and d

θ/dζ

(partial derivative symbols are obviously no longer needed). Thus:

∂T

∂t

dθ

dζ

−ζ

(

−T

)

dθ

dζ

(3.1.6)

132 Energy transfer processes in planetary bodies

Box 3.1

Continued

(you should work out the algebra to see where this comes from). Similarly:

∂T

∂x

dθ

dζ

−T

√

κt

dθ

dζ

(3.1.7)

and:

∂

∂x

dζ





dζ

−T

4κt

dζ

. (3.1.8)

We now substitute (3.1.6) and (3.1.8) into (3.1.1) and, after simplifying, we get:

−2ζ

dθ

dζ

∂

∂ζ

. (3.1.9)

In order to solve this ODE we introduce the auxiliary variable φ:

φ =

dθ

dζ

(3.1.10)

and its ﬁrst derivative:

dφ

dζ

. (3.1.11)

Substituting into (3.1.9) and rearranging:

dφ

=−2ζ dζ . (3.1.12)

In case you lost track, we have converted the original PDE, equation (3.1.1), into the very simple ordinary

differential equation (3.1.12). It may be a good idea to go back and review all of the steps involved in this

rather clever trick. The solution to (3.1.12) is, simply:

φ = ce

−ζ

, (3.1.13)

where c is an integration constant. We can now do away with the auxiliary variable φ, and we get:

dθ

dζ

=ce

−ζ

. (3.1.14)

In order to ﬁnd the value of the constant c we integrate between the two conditions speciﬁed by equations

(3.1.5):



θ=1

θ=0

dθ = c



ζ =∞

ζ =0

−ζ

dζ . (3.1.15)

133 3.2 Heat transport by diffusion

Box 3.1

Continued

The integral on the right-hand side of this equation appears in a function called the error function (see, for

example, Weisstein, 2003, p. 933), symbolized erf, and deﬁned by:

erf

(

)

√



η=ζ

η=0

−η

dη, (3.1.16)

where η is a dummy variable of integration. Equation (3.1.15) then becomes:

1 =

√

erf

(

∞

)

. (3.1.17)

The values of the error function are tabulated (see Beyer, 1987) and also available in Maple. We ﬁnd that erf

(∞) =1, so that the integration constant is:

c =

√

. (3.1.18)

Integrating (3.1.14) with this constant:



dθ =

√



−η

dη (3.1.19)

or:

θ = erf

(

)

. (3.1.20)

This is the solution that we seek: the function that describes temperature distribution in a semi-inﬁnite half

space cooled by diffusion and without a heat source. It is straightforward to substitute the deﬁnitions of the

non-dimensional variables ζ and θ (equations (3.1.3) and (3.1.4)) into equation (3.1.20) and obtain the

equation for T as a function of x and t:

T =

(

−T

)

erf



√

κt



. (3.1.21)

Solution of the heat diffusion equation in a semi-inﬁnite half space is a relatively simple exercise but it is

not necessarily the correct answer to every problem in heat diffusion. The classic text by Carlsaw and Jaeger

(1959) develops a large number of solutions of the heat equation with different geometries and boundary

conditions, many of them applicable to problems in the Earth and planetary sciences.

3.2.3 The physical meaning of diffusivity

Diffusivity measures the efficiency with which interactions at the microscopic level are

able to eliminate the potential gradient that drives the flow. In the expression κ = k/(ρc

)

we can think of the product ρc

as the thermal inertia of the material. The larger the

value of this product, the more difficult it is to change the thermal state of the material,

thus slowing down heat diffusion. All diffusivities (for heat, momentum and mass) have

dimensions of [L]

[T ]

−1

, as should be evident from the fact that they have to satisfy the

134 Energy transfer processes in planetary bodies

0 0.2

Increasing distance

from cooled boundary

Increasing time

0.4 0.6

0.8 1

cooling

x=2 t

x= t

x=4 t

Thermal

boundary

layer

Fig. 3.5

The error function θ = erf (ζ ). Note that the axes for the independent and dependent variables are switched relative

to the mathematical convention, as is common practice in the planetary sciences when representing temperature

as a function of depth. The similarity variable ζ =x/2

√

κt tracks the propagation of a temperature perturbation in

space and time, as shown by the arrows. The magnitude of the perturbation is given by the distance from the

non-dimensional temperature θ = 1. The thermal boundary layer can be (arbitrarily) deﬁned at ζ = 1, for

which θ ≈0.84.

same relationship between spatial and temporal derivatives given by equation (3.15). Other

physical quantities that appeared in the derivation of the diffusion equation, most notably

temperature and energy, do not show up in the dimension of diffusivity. The dimension of

diffusivity tells us that in all diffusive processes length scales with the square root of time.

Consider the simple case of diffusion in one dimension (i.e. ∂T /∂y =∂T /∂z =0) without

a source (α =0). In Box 3.1 we solve this equation for a semi-infinite half space. Let us look

at the general solution in terms of the non-dimensional similarity variable ζ = x/2

√

(κt)

and non-dimensional temperature, θ. The error function, which is plotted in terms of these

variables in Fig. 3.5, describes the propagation of a temperature perturbation (θ =0) imposed

at the boundary of a half space, for which the uniform initial condition is θ =1. The distance

between the curve and the right vertical axis (θ =1) is the amount of cooling, which increases

as ζ decreases. If we fix time, then a displacement towards ζ =0 represents a displacement

towards the boundary of the half space (x = 0), which, for finite time, is the only point for

which θ = 0. Alternatively, if we consider a fixed point in space, then time, and with it the

amount of cooling, increases towards ζ = 0, and cooling to θ = 0 for an arbitrary value of

x = 0 can be achieved only after an infinite amount of time (of course, this is true for a

“semi-infinite” half space, as complete cooling of a finite body such as a planet, a dike or

lava flow will of course always be accomplished in a finite interval of time).

135 3.2 Heat transport by diffusion

We can now define a value of ζ that gives a scale estimate for the propagation of a

temperature perturbation. As an example, if we set x =

√

(κt), we get ζ =

and θ ≈ 0.52

(Fig. 3.5). Temperature at this point has changed substantially from its initial value, but the

effect of the temperature perturbation extends considerably more than this. If we choose

ζ =2, we find that θ =erf (2) =0.995, which means that the temperature perturbation at this

point is insignificant. An intermediate value between these two would give a good estimate

of how far the temperature perturbation has propagated. The choice is always arbitrary, and

in this book we will choose ζ = 1, for which θ = erf (1) = 0.843 (Fig. 3.5). Although the

temperature perturbation can still be seen for ζ>1, most of the cooling has taken place in

the interval 0 <ζ <1. This gives us the following relationship between the characteristic

diffusion length, λ, and characteristic diffusion time, τ :

λ ∼2

√

κτ, (3.16)

where the symbol “∼” means “of the same order of magnitude as”. Given that θ is non-

dimensional temperature, this relationship holds regardless of the absolute magnitude of

the initial temperature perturbation.

The physical meaning of equation (3.16) is that a temperature perturbation of any size

will travel a distance of order λ in a time τ . The numerical value of the diffusivity provides

a scale estimate for heat conduction. Consider for example the oceanic lithosphere. It forms

by diffusive cooling of mantle asthenosphere as it moves away from spreading centers. We

can therefore take the age of the ocean floor as the cooling time of the lithosphere. Using an

average age of the ocean floor of 100 million years, and a typical value of heat diffusivity

for rocks of 10

−6

−1

, we calculate from (3.16) an average thickness for the oceanic

lithosphere of ∼112 km, which is certainly of the right order of magnitude. Similarly, we

can estimate that the characteristic cooling time of a 10 m-thick ash flow (assuming that we

can neglect circulation of hydrothermal fluids) is ∼10 months, which means that, even if

the surface of the ash flow cools to room temperature much sooner than this, we can expect

that temperatures deep inside the ash flow will remain noticeably above room temperature

for about a year. Recall that in diffusive processes time scales with the square of length.

Thus, a 1-m thick lava flow cools 100 times faster than a 10-m thick lava flow – it will have

cooled in ∼3 days.

These numerical examples illustrate a powerful method for obtaining approximate

answers called scale analysis or scaling. The goal of scaling is to obtain a numerical value

that we can be reasonably certain is within an order of magnitude of the exact answer that

we seek. In many cases obtaining the exact answer may be computationally very intensive

and a scale analysis will return most of the information that we seek for a minimum invest-

ment of time and effort. This is especially true in planetary sciences, where the complexity

of the processes involved and the large uncertainties in the values of material properties

and intensive variables may not always render an “exact” calculation any more informative

than a scale analysis. The spirit of scale analysis is to capture the essential physics of the

phenomenon that we are studying, even as many of the details may remain uncertain or

unknown. One must be careful not to miss some important information, though. In the words

of Albert Einstein, “Everything should be made as simple as possible, but not simpler”.

3.2.4 The diffusive thermal boundary layer

The thickness of the half space that has been affected by the temperature perturbation (as

defined by equation (3.16)) is called the thermal boundary layer. Most of the temperature

136 Energy transfer processes in planetary bodies

difference between the interior of the half space and the boundary occurs across the thermal

boundary layer (Fig. 3.5). A physical interpretation of the diffusive thermal boundary layer

is that the material beyond it does not “know” that it will be changing its temperature and,

therefore, does not contribute internal energy to the heat flow across the boundary of the

half space. Heat flow across the boundary comes only from cooling of material inside the

thermal boundary layer.

We seek an expression for the rate of heat loss across the thermal boundary layer. We

begin from the solution to the heat diffusion equation in a semi-infinite half space (equation

(3.1.21)):

T =

(

−T

)

erf



√

κt



(3.17)

and using the similarity variable ζ defined by equation (3.1.3) we get, from equation 3.17

and the chain rule:

(

−T

)



dζ

erf

(

)



dζ

−T

√

κt

dζ

erf

(

)

. (3.18)

From the definition of the error function, equation (3.1.16), it follows that:

dζ

erf

(

)

√

−ζ

(3.19)

so that the temperature gradient is given by:

−T

√

πκt

−ζ

. (3.20)

Since at the boundary ζ =x= 0, the heat flux across the boundary of the half space is

given by:

=−k



=−k

−T

√

πκt

. (3.21)

When applied to surface heat flux in planetary bodies, the negative sign in this equation can

be the source of some consternation. If we measure depths in a planet as positive quantities,

then heat flux at the planet’s surface is a negative quantity, as this equation, which is just

Fourier’s law, shows. It is customary, however, to express planetary heat flux as a positive

quantity, which is what we have been doing throughout this book. In order to follow this

near universal practice, which we will continue to do, it is necessary to drop the negative

sign from equation (3.21), which thus becomes:

−T

√

πκt

(3.22)

and it is also necessary to remember that we are inverting the sign of the planetary heat

flux, and adjust the sign in other equations, as needed. As the boundary layer cools and

thickens the thermally unperturbed material recedes progressively farther away (Fig. 3.5).

By Fourier’s law (equation (3.5)) heat flux across the boundary of the half space must also

decrease with time, as shown by equation (3.22).

137 3.3 Heat diffusion and cooling of planetary bodies

3.3 Heat diffusion and cooling of planetary bodies

3.3.1 Lord Kelvin and the thermal structure of the Earth

Students of the Earth and planetary sciences are familiar with the nineteenth-century con-

troversy about the age of the Earth between Lord Kelvin, who calculated an age of the

order of 10

–10

years, and his contemporary geologists, who insisted on much older ages.

Insofar as some geologists advocated an undetermined, perhaps infinitely old age, and by

implication a steady-state Earth, this was an argument about the Second Law of Thermo-

dynamics. Kelvin, who was one of the foundational figures of classical thermodynamics,

correctly pointed out that a steady state Earth would violate the Second Law (Chapter 4).

I believe that the outcome of this argument was split. Kelvin’s was a quantitative estimate

based on a rigorous application of physics, and in this sense it was better science than

anything that nineteenth-century geology could offer. Lord Kelvin was also correct that the

Earth is not infinitely old. The age of the planet, however, is one to two orders of magnitude

greater than suggested by his calculations. As it turns out, Lord Kelvin’s error provides a

crucial clue about the interior structure of the Earth, which was not fully recognized until

the 1960s (although some scientists had grasped the significance of his error much earlier;

see England et al., 2007, for an account of one such case). Our focus here is to see what

we can learn about the interior of the Earth and other terrestrial planets from the physics of

heat diffusion.

Let us assume that we are ignorant of the interior structure of the Earth and that we

assume that it is a largely solid sphere in which heat moves by diffusion. Because there is

measurable heat flow at the surface, we can also conclude that, in the absence of an active

heat source, the Earth is cooling down from an initially hotter state. This was, more or less,

Lord Kelvin’s view of the Earth. Let us now add a crucial piece of information that was

unknown to him: the correct age of the Earth, 4.56 Ga ≈1.44 ×10

s. From equation (3.16),

with κ =10

−6

−1

, we can conclude that only the outermost 750 km or so of the Earth

have cooled down. Deeper than that the Earth should still preserve temperatures comparable

to its primordial temperature, because there has not been enough time for diffusion to extract

heat from those depths. The observed heat flow at the surface of the planet must therefore

derive exclusively from cooling of the outermost 750-km thick layer. If we know the age of

the Earth then we can use the diffusion equation to do one of two things: assume an initial

temperature for the Earth and calculate the surface heat flux, or use the measured heat flux

and calculate the initial temperature of the Earth. We can also do what Lord Kelvin did:

use measured heat flux and an estimate of the Earth’s initial temperature to calculate the

length of time for which the Earth has been cooling down. We will now perform all three

calculations and show that they all lead to erroneous results. More accurately, they lead to

three ways of looking at the same result, namely, that heat diffusion by itself cannot explain

the thermal structure of the Earth.

We model the diffusively cooled Earth as a semi-infinite half space. This is acceptable

because, over the age of the Earth, a layer with thickness of only about one tenth of the

Earth’s radius can cool by diffusion, so that most of the planet effectively lies infinitely

distant from the surface in terms of heat transfer. We set x = 0 at the surface of the Earth

and x>0 towards the center of the Earth. We assume that the Earth is initially at uniform

temperature T

, and that temperature at the surface has a constant value T

for all times.