# Heat Conduction Cheat Sheet

I am dumping some equations here I need now and then! The sections about 3-dimensional temperature waves summarize what is described at length in the second part of this post.

Temperature waves are interesting for simulating yearly and daily oscillations in the temperature below the surface of the earth or near wall/floor of our ice/water tank. Stationary solutions are interesting to assess the heat transport between pipes and the medium they are immersed in, like the tubes making up the heat exchanger in the tank or the solar/air collector.

~

## Heat equation – conservation of energy [Top]

Energy is conserved locally. It cannot be destroyed or created, but it it is also not possible to remove energy in one spot and make it reappear in a different spot. The energy density η in a volume element can only change because energy flows out of this volume, at a flow density j (energy per area and time).

$\frac{\partial \eta}{\partial t} + \frac{\partial \vec{j}}{\partial\vec{r}} = 0$

In case of heat energy, the sensible heat energy ‘contained’ in a volume element is the volume times mass density ρ [kg/m3] times specific heat c [J/kgK] times the temperature difference in K (from a ‘zero point’). The flow of heat energy is proportional to the temperature gradient (with constant λ – heat conductivity [J/mK], and heat flows from hot to colder spots.

$\rho c \frac{\partial T}{\partial t} + \frac{\partial}{\partial\vec{r}} (- \lambda \frac{\partial T}{\partial\vec{r}}) = 0$

Re-arranging and assuming that the three properties ρ, c, and λ are constant in space and time, they can be combined into a single property called thermal diffusivity D

$D = \frac{\lambda}{\rho c}$

$\frac{\partial T}{\partial t} = D \frac{\partial}{\partial\vec{r}} \frac{\partial T}{\partial\vec{r}} = D \Delta T$

In one dimensions – e.g. heat conduction to/from an infinite plane –  the equation is

$\frac{\partial T}{\partial t} = D \frac{d^{2} T}{d x^{2}}$

~

## 1D solution – temperature waves in one dimension [Top]

I covered it already here in detail. I’m using complex solutions as some manipulations are easier to do with the exponential functions than with trigonometric functions, keeping in mind we are finally interested in the real part.

Boundary condition – oscillating temperature at the surface; e.g. surface temperature of the earth in a year. Angular frequency ω is 2π over period T (e.g.: one year)

$T(t,0) = T_0 e^{i \omega t}$

Ansatz: Temperature wave, temperature oscillating with ω in time and with to-be-determined complex β in space.

$T(t,x) = T_0 e^{i \omega t + \beta x}$

Plugging into 1D heat equation, you get β as a function of ω and the properties of the material:

$i \omega = D \beta^2$
$\beta = \pm \sqrt{\frac{i \omega}{D}} = \pm \sqrt{i} \sqrt{\frac{\omega}{D}} = \pm (1 + i){\sqrt 2} \sqrt{\frac{\omega}{D}} = \pm (1 + i) \sqrt{\frac{\omega}{2D}}$

The temperature should better decay with increasing x – only the solution with a negative sense makes sense, then $T(\infty) = T_0$. The temperature well below the surface, e.g. deep in the earth, is the same as the yearly average of the air temperature (neglecting the true geothermal energy and related energy flow and linear temperature gradient).

Solution – temperature as function of space and time:

$T(t,x) = T_0 e^{i \omega t - (1 + i) \sqrt{\frac{\omega}{2D}} x} = T_0 e^{i (\omega t - \sqrt{\frac{\omega}{2D}} x)} e^{-\sqrt{\frac{\omega}{2D}} x}$

Introducing parameter k:

$\sqrt{\frac{\omega}{2D}} = k$

Concise version of the solution function:

$T(t,x) = T_0 e^{i (\omega t - kx)} e^{-kx}$

Strip off the real part:

$Re(T(t,x)) = T_0 cos(\omega t - kx) e^{-kx}$

Relations connecting the important wave parameters:

$\tau = \frac {2 \pi}{\omega}$
$\lambda = \frac {2 \pi}{k}$

~

## ‘Helpers’ for the 3D case (spherical) [Top]

Basic stuff

$r = \sqrt{x^2 + y^2 + z^2}$
$\frac{\partial r}{\partial \vec{r}} = (\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z})\sqrt{x^2 + y^2 + z^2} = \frac{\vec{r}}{r}$
$\frac{\partial \vec{r}}{\partial \vec{r}} = (\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial z}{\partial z})(x,y,z) = 3$
$\Delta T = (\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2})T(r)$

Inserting, to obtain a nicely looking Laplacian in spherical symmetry

$\Delta T = \frac{\partial}{\partial\vec{r}} \frac{\partial}{\partial\vec{r}} T(\sqrt{x^2 + y^2 + z^2}) = \frac{\partial}{\partial\vec{r}} \frac{\partial r}{\partial\vec{r}} (\frac{dT}{dr}) = \frac{\partial}{\partial\vec{r}} (\frac{\vec{r}}{r} \frac{dT}{dr})$
$= \frac{3}{r} \frac{dT}{dr} - \frac{1}{r^2} \frac{\partial r}{\partial\vec{r}} \vec{r} \frac{dT}{dr} + \frac{\vec{r}}{r} \frac{\vec{r}}{r} \frac{d^2 T}{dr^2}$
$= \frac{3}{r} \frac{dT}{dr} - \frac{1}{r} \frac{dT}{dr}+ \frac{d^2 T}{dr^2} = \frac{2}{r} \frac{dT}{dr} + \frac{d^2 T}{dr^2}$
$= \frac{1}{r}(\frac{dT}{dr} + \frac{dT}{dr} + r \frac{d^2T}{dr^2}) = \frac{1}{r} \frac{d}{dr} (T + r \frac{dT}{dr}) = \frac{1}{r} \frac{d^2}{dr^2}(rT)$

~

## ‘Helpers’ for the 2D case (cylindrical) [Top]

Basic stuff

$r = \sqrt{x^2 + y^2}$
$\frac{\partial r}{\partial \vec{r}} = (\frac{\partial}{\partial x}$

~

## 3D stationary solution- sphere [Top]

The time derivative is zero, so the Laplacian is zero:

$\frac{\partial T}{\partial t} = 0 = \Delta T(t, r) = \frac{1}{r} \frac{d^2}{dr^2}(rT)$

Ansatz, guessing something simple

$T(r) = \frac{A + Br}{r} = \frac{A}{r} + B$

Boundary conditions, as for the 1D case:

$T(R_0) = T_0$
$T(R_E) = T_E$

Plugging in – getting functions for all r:

$T(r) = \frac{1}{R_0 - R_E} [R_E T_E(\frac{R_0}{r} - 1) + R_0 T_0 (1 - \frac{R_E}{r}]$

$|j(r)| = \lambda \frac{1}{R_0 - R_E} \frac{1}{r^2} [R_E T_E R_0 - R_0 T_0 R_E ]$

At the surface:

$|j(R_0)| = \lambda \frac{1}{R_0 - R_E} \frac{R_E}{R_0} [T_E - T_0 ]$

~

## 2D stationary solution – cylinder, pipe [Top]

Cylindrical Laplacian is zero

$\frac{1}{r} \frac{d}{dr} (r \frac{dT}{dr}) = 0$

Same boundary conditions, plugging in

$r \frac{dT}{dr} = A$
$dT = A \frac {dr}{r}$

$\int_{T}^{T_0} dT = A \int_{R_0}^{r} \frac {dr}{r}$
$T(r) = T_0 + A \ln{(\frac{r}{R_0})} = T_0 + A (\ln{r} - \ln{R_0})$
$T(R_E) = T_E = T_0 + A \ln{(\frac{R_E}{R_0})}$
$A = \frac{T_E - T_0}{\ln{(\frac{R_E}{R_0})}}$

Solutions for temperature and energy flow at any r:

$T(r) = T_0 + (T_E - T_0) \frac{\ln{(\frac{r}{R_0})}}{\ln{(\frac{R_E}{R_0})}}$
$|\vec{j(r)}| = |\frac {1}{r} \lambda \frac{T_E - T_0}{\ln{(\frac{R_E}{R_0})}}|$

Expressing r in terms of distance from the surface, $\Delta r = r - R_0$

$|\vec{j(r)}| = |\frac {1}{\Delta r + R_0} \lambda \frac{T_E - T_0}{\ln{(\frac{R_1}{R_0})}}|$

~

## Comparison of overall heat flow: 1D versus 2D [Top]

j is the energy flow per area, and the area traversed by the flow depends on geometry. in the 1D case the area is always the same area, equal to the area of the plane. For a cylinder, the area increases with r.

The integrated energy flow J for a plate with area F is

$J_{Plate} = F \lambda \frac{|T_E - T_0|}{\Delta x}$

If the two temperatures are given, J decreases linearly with increasing thickness of the cylindrical ‘shell’, e.g. a growing layer of ice.

For a cylinder of length l the energy flow J is…

$J_{Cyl} = 2 \pi l r |\frac {1}{r} \lambda \frac{T_E - T_0}{\ln{(\frac{R_E}{R_0})}}|$
$= 2 \pi l \lambda |\frac{T_E - T_0}{\ln{(\frac{R_E}{R_0})}}|$ \par

Factor r has been cancelled, and the for given temperatures J is only decreasing linearly with increasing outer radius $R_E$. That’s why vendors of plate heat exchangers (in vessels with phase change material) worry more about a growing layer of sold material than user for e.g. ‘ice on coil’ I quoted a related research paper on ‘ice storage powered’ heat pump system in this post – they make exactly this point and provide some data. In addition to conduction also convection at both sides of the heat exchanger should be taken into account, too, in a ‘serial connection’ of heat transferring components.

# Rowboats, Laser Pulses, and Heat Energy (Boring Title: Dimensional Analysis)

Dimensional analysis means to understand the essentials of a phenomenon in physics and to calculate characteristic numbers – without solving the underlying, often complex, differential equation. The theory of fluid dynamics is full of interesting dimensionless numbers –  Reynolds Number is perhaps most famous.

In the previous post on temperature waves I solved the Heat Equation for a very simple case, in order to answer the question How far does solar energy get into ground in a year? Reason: I have been working on simulations of our heat pump system since a few years. This also involves heat transport between the water/ice tank and ground. If you set out to simulate a complex phenomenon you have to make lots of assumptions about materials’ parameters, and you have to simplify the system and equations you use for modelling the real world. You need a way of cross-checking if your results sound plausible in terms of orders of magnitude. So my goal has been to find yet another method to confirm assumptions I have made about the thermal properties of ground elsewhere.

Before I am going to revisit heat transport, I’ll try to explain what dimensional analysis is – using the best example I’ve ever seen. I borrow it from theoretical physicist – and awesome lecturer – David Tong:

How does the speed of a rowing boat depend in the number of rowers?

References: Tong’s lecture called Dynamics and Relativity (Chapter 3), This is the original paper from 1971 Tong quotes: Rowing: A similarity analysis.

The boat experiences a force of friction in water. As for a car impeded by the friction of the surrounding air, the force of friction depends on velocity.

Force is the change of momentum, momentum is proportional to mass times velocity. Every small ‘parcel’ of water carries a momentum proportional to speed – so force should at least be proportional to one factor of v. But these parcel move at a speed v, so the faster they move the more momentum is exchanged with the boat; so there has to be a second factor of v, and force is proportional to the square of the speed of the boat.

The larger the cross-section of the submerged part of the boat, A, the higher is the number of collisions between parcels of water and the boat, so putting it together:

$F \sim v^{2}A$

Rowers need to put in power to compensate for friction. Power is energy per time, and Energy is force times distance. Since distance over time is velocity, thus power is also force times velocity.

So there is one more factor of v to be included in power:

$P \sim v^{3}A$

For the same reason wind power harvested by wind turbines is proportional to the third power of wind speed.

A boat does not sink because downward gravity and upward buoyancy just compensate each other; buoyancy is the weight of the volume of water displaced. The heavier the load, the more water needs to be displaced. The submerged volume of the boat V is proportional to the weight of the rowers, and thus to their number N if the mass of the boat itself is negligible:

$V \sim N$

The volume of something scales with the third power of its linear dimensions – think of a cube or a sphere; so the surface area scales with the square of the length, and the cross-section A scales with V – and thus with N:

$A \sim N^{\frac{2}{3}}$

Each rower contributes the same share to the total rowing power, so:

$P \sim N$

Inserting for A in the first expression for P:

$P \sim v^{3} N^{\frac{2}{3}}$

Eliminating P as it has been shown to be proportional to N:

$N \sim v^{3} N^{\frac{2}{3}}$
$v^{3} \sim N^{\frac{1}{3}}$
$v \sim N^{\frac{1}{9}}$

… which is in good agreement with measurements according to Tong.

Heat Transport and Characteristic Lengths

In the last post I’ve calculated characteristic lengths, describing how heat is slowly dissipated in ground: 1) The wavelength of the damped oscillation and 2) the run-out length of the enveloping exponential function.

Both are proportional to the square root of a simple number:

$l \sim \sqrt{D \tau}$

… the factor of proportionality being ‘small’ on a logarithmic scale, like π or 2 or their inverse. τ is the period, and D was a number expressing how well the material carries away heat energy.

There is another ‘simple’ scenario that also results in a length scale described by
$\sqrt{D \tau}$ times a small number: If you deposit a confined ‘lump of heat’, a ‘point heat’ it will peter out and the average width of the lump after some time τ is about this length as well.

Using very short laser pulse to heat solid material is very close to depositing ‘point heat’. Decades ago I worked with pulsed excimer lasers, used for ablation (‘shooting off) material from ceramic targets.This type of lasers is used in eye surgery today:

Heat is deposited in nanosecond pulses, and the run-out length of the heat peak in the material is about $\sqrt{D \tau}$ with tau being equal to the very short laser’s pulse length of several nanoseconds. As the pulse duration is short, the penetration depth is short as well, and tissue is ‘cut’ precisely without heating much of the underlying material.

So this type of $\sqrt{D \tau}$ length is not just a result of a calculation for a specific scenario, but it rather seems to encompass important characteristics of heat conduction as such.

The unit of D is area over time, m2/s. If you accept the heat equation as a starting point, analysing the dimensions involved by counting x and t you see that D has to contain two powers of x and one of t. Half of applied physics and engineering is about getting units right.

But I pretend I don’t even know the heat equation and ‘visualize’ heat transport in this way: ‘Something’ – like heat energy – is concentrated in space and closely peters out. The spreading out is faster, the more concentrated it is. A thin needle-like peak quickly becomes a rounded hill, and then is flattened gradually. Concentration in space means curvature. The smaller the space occupied by the lump of heat is, the smaller its radius, the higher its curvature as curvature is the inverse of the radius of a tangential circular path.

I want to relate curvature to the change with time. Change in time has to be measured in units including the inverse of time, curvature is the inverse of space. Equating those, you have to come with something including the square of spatial dimension and one temporal dimension – something like D [m2/s].

How to get a characteristic length from this? D has to be multiplied by a characteristic time, and then we need to take a the square root. So we need to put in some characteristic time, that’s a property of the specific system investigated and not of the equation – like the yearly period or the laser pulse. And the resulting length is exactly that $l \sim \sqrt{D \tau}$ that shows up in any of of the solutions for specific scenarios.

_________________________________

The characteristic width of the spreading lump of heat is visible in the so-called Green’s functions. These functions described a system’s response to a ‘source’ which resemble a needle-like peak in time. In this case it is a Gaussian functions with a ‘width’ $\sim \sqrt{D \tau}$. See e.g. equation (14) on PDF-page 14 of these lecture notes.

Penetration depth of excimer lasers in human tissue – in this book the square root D times tau formula is used and depths are calculated to be equal to several 10 micrometers.

# Temperature Waves and Geothermal Energy

Nearly all of renewable energy exploited today is, in a sense, solar energy. Photovoltaic cells convert solar radiation into electricity, solar thermal collectors heat hot water. Plants need solar power for photosynthesis, for ‘creating biomass’. The motion of water and air is influenced by the fictitious forces caused by the earth’s rotation, but by temperature gradients imposed by the distribution of solar energy as well.

Also geothermal heat pumps with ground loops near the surface actually use solar energy deposited in summer and stored for winter – that’s why I think that ‘geothermal heat pumps’ is a bit of a misnomer.

Collector (heat exchanger) for brine-water heat pumps.

Within the first ~10 meters below the surface, temperature fluctuates throughout the year; at 10m the temperature remains about constant and equal to 10-15°C for the whole year.

Only at higher depths the flow of ‘real’ geothermal energy can be spotted: In the top layer of the earth’s crust the temperatures rises about linearly, at about 3°C (3K) per 100m. The details depend on geological peculiarities, it can be higher in active regions. This is the energy utilized by geothermal power plants delivering electricity and/or heat.

Geothermal gradient adapted from Boehler, R. (1996). Melting temperature of the Earth’s mantle and core: Earth’s thermal structure. Annual Review of Earth and Planetary Sciences, 24(1), 15–40. (Wikimedia, user Bkilli1). Geothermal power plants use boreholes a few kilometers deep.

This geothermal energy originates from radioactive decays and from the violent past of the primordial earth: when the kinetic energy of celestial objects colliding with each other turned into heat.

The flow of geothermal energy per area directed to the surface, associated with this gradient is about 65 mW/m2 on continents:

Global map of the flow of heat, in mW/m2, from Earth’s interior to the surface. Davies, J. H., & Davies, D. R. (2010). Earth’s surface heat flux. Solid Earth, 1(1), 5-24. (Wikimedia user Bkilli1)

Some comparisons:

• It is small compared to the energy from the sun: In middle Europe, the sun provides about 1.000 kWh per m2 and year, thus 1.000.000Wh / 8.760h = 144W/m2 on average.
• It also much lower than the rule-of-thumb power of ‘flat’ ground loop collectors – about 20W/m2
• The total ‘cooling power’ of the earth is several 1010kW: Would the energy not be replenished by radioactive decay, the earth would lose a some seemingly impressive 1014kWh per year, yet this would result only in a temperature difference of ~10-7°C (This is just a back-of-the-envelope check of orders of magnitude, based on earth’s mass and surface area, see links at the bottom for detailed values).

The constant energy in 10m depth – the ‘neutral zone’ – is about the same as the average temperature of the earth (averaged over one year over the surface of the earth): About 14°C. I will show below that this is not a coincidence: The temperature right below the fluctuating temperature wave ‘driven’ by the sun has to be equal to the average value at the surface. It is misleading to attribute the 10°C in 10m depths to the ‘hot inner earth’ only.

In this post I am toying with theoretical calculations, but in order not so scare readers off too much I show the figures first, and add the derivation as an appendix. My goal is to compare these results with our measurements, to cross-check assumptions for the thermal properties of ground I use in numerical simulations of our heat pump system (which I need for modeling e.g. the expected maximum volume of ice)

1. The surface temperature varies periodically in a year, and I use maximum, minimum and average temperature from our measurements, (corrected a bit for the mild last seasons). These are daily averages as I am not interested in the daily temperature changes between and night.
2. A constant geothermal flow of 65 mW/m2 is superimposed to that.
3. The slow transport of solar energy into ground is governed by a thermal property of ground, called the thermal diffusivity. It describes ‘how quickly’ a lump of heat deposited will spread; its unit is area per time. I use an assumption for this number based on values for soil in the literature.

I am determining the temperature as a function of depth and of time by solving the differential equation that governs heat conduction. This equation tells us how a spatial distribution of heat energy or ‘temperature field’ will slowly evolve with time, given the temperature at the boundary of the interesting part of space in question – in this case the surface of the earth. Fortunately, the yearly oscillation of air temperature is about the simplest boundary condition one could have, so you can calculate the solution analytically.
Another nice feature of the underlying equation is that it allows for adding different solutions: I can just add the effect of the real geothermal flow of energy to the fluctuations caused by solar energy.

The result is a  ‘damped temperature wave’; the temperature varies periodically with time and space: The spatial maximum of temperature moves from the surface to a point below and back: In summer (beginning of August) the measured temperature is maximum at the surface, but in autumn the maximum is found some meters below – heat flows back from ground to the surface then:

Calculated ground temperature, based on measurements of the yearly variation of the temperature at the surface and an assumption of the thermal properties of ground. Calculated for typical middle European maximum and minimum temperatures.

This figure is in line with the images shown in every textbook of geothermal energy. Since the wave is symmetrical about the yearly average, the temperature in about 10m depth, when the wave has ‘run out’, has to be equal to the yearly average at the surface. The wave does not have much chance to oscillate as it is damped down in the middle of the first period, so the length of decay is much shorter than the wavelength.

The geothermal flow just adds a small distortion, an asymmetry of the ‘wave’. It is seen only when switching to a larger scale.

Some data as in previous plot, just extended to greater depths. The geothermal gradient is about 3°C/100m, the detailed value being calculated from the value of thermal conductivity also used to model the fluctuations.

Now varying time instead of space: The higher the depth, the more time it takes for ground to reach maximum temperature. The lag of the maximum temperature is proportional to depth: For 1m difference in depth it is less than a month.

Temporal change of ground temperature at different depths. The wave is damped, but other simply ‘moving into the earth’ at a constant speed.

Measuring the time difference between the maxima for different depths lets us determine the ‘speed of propagation’ of this wave – its wavelength divided by its period. Actually, the speed depends in a simple way on the thermal diffusivity and the period as I show below.

But this gives me an opportunity to cross-check my assumption for diffusivity: I  need to compare the calculations with the experimentally determined delay of the maximum. We measure ground temperature at different depths, below our ice/water tank but also in undisturbed ground:

Temperature measured with Pt1000 sensors – comparing ground temperature at different depths, and the related ‘lag’. Indicated by vertical dotted lines, the approximate positions of maxima and minima. The lag is about 10-15 days.

The lag derived from the figure is in the same order as the lag derived from the calculation and thus in accordance with my assumed thermal diffusivity: In 70cm depth, the temperature peak is delayed by about two weeks.

___________________________________________________

Appendix: Calculations and background.

I am trying to give an outline of my solution, plus some ‘motivation’ of where the differential equation comes from.

Heat transfer is governed by the same type of equation that describes also the diffusion of gas molecules or similar phenomena. Something lumped together in space slowly peters out, spatial irregularities are flattened. Or: The temporal change – the first derivative with respect to time – is ‘driven’ by a spatial curvature, the second derivative with respect to space.

$\frac{\partial T}{\partial t} = D\frac{\partial^{2} T}{\partial x^{2}}$

This is the heat transfer equation for a region of space that does not have any sources or sinks of heat – places where heat energy would be created from ‘nothing’ or vanish – like an underground nuclear reaction (or freezing of ice). All we know about the material is covered by the constant D, called thermal diffusivity.

The equation is based on local conservation of energy: The energy stored in a small volume of space can only change if something is created or removed within that volume (‘sources’) or if it flows out of the volume through its surface. This is a very general principles applicable to almost anything in physics. Without sources or sinks, this translates to:

$\frac{\partial [energy\,density]}{\partial t} = -\frac{\partial \overrightarrow{[energy\,flow]}}{\partial x}$

The energy density [J/m3] stored in a volume of material by heating it up from some start temperature is proportional to temperature, proportionality factors being the mass density ρ [kg/m3] and the specific heat cp [J/kg] of this material. The energy flow per area [W/m2] is typically nearly proportional to the temperature gradient, the constant being heat conductivity κ [W/mK]. The gradient is the first-order derivative in space, so inserting all this we end with the second derivative in space.

All three characteristic constants of the heat conducting material can be combined into one – the diffusivity mentioned before:

$D = \frac{\kappa }{\varrho \, c_{p} }$

So changes in more than one of these parameters can compensate for each other; for example low density can compensate for low conductivity. I hinted at this when writing about heat conduction in our gigantic ice cube: Ice has a higher conductivity and a lower specific heat than water, thus a much higher diffusivity.

I am considering a vast area of ground irradiated by the sun, so heat conduction will be one-dimensional and temperature changes only along the axis perpendicular to the surface. At the surface the temperature varies periodically throughout the year. t=0 is to be associated with beginning of August – our experimentally determined maximum – and the minimum is observed at the beginning of February.

This assumption is just the boundary condition needed to solve this partial differential equation. The real ‘wavy’  variation of temperature is closed to a sine wave, which makes the calculation also very easy. As a physicist I have trained to used a complex exponential function rather than sine or cosine, keeping in mind that only real part describes the real world. This a legitimate choice, thanks to the linearity of the differential equation:

$T(t,x=0) = T_{0} e^{i\omega t}$

with ω being the angular frequency corresponding to one year (2π/ω = 1 year).

It oscillates about 0, with an amplitude of T0. But after all, the definition of 0°C is arbitrary and – again thanks to linearity – we can use this solution and just add a constant function to shift it to the desired value. A constant does neither change with space or time and thus solves the equation trivially.

If you have more complicated sources or sinks, you would represent those mathematically as a composition of simpler ‘sources’, for each of which you can find a quick solution and then add up add the solutions, again thanks to linearity. We are lucky that our boundary condition consist just of one such simple harmonic wave, and we guess at the solution for all of space, adding a spatial wave to the temporal one.

So this is the ansatz – an educated guess for the function that we hope to solve the differential equation:

$T(t,x) = T_{0} e^{i\omega t + \beta x}$

It’s the temperature at the surface, multiplied by an exponential function. x is positive and increasing with depth. β is some number we don’t know yet. For x=0 it’s equal to the boundary temperature. Would it be a real, negative number, temperature would decrease exponentially with depth.

The ansatz is inserted into the heat equation, and every differentiation with respect to either space or time just yields a factor; then the exponential function can be cancelled from the heat transfer equation. We end up with a constraint for the factor β:

$i\omega = D\beta^{2}$

Taking the square root of the complex number, there would be two solutions:

$\beta=\pm \sqrt{\frac{\omega}{2D}}(1+i))$

β has a real and an imaginary part: Using it in T(x,t) the real part corresponds to exponential ‘decay’ while the imaginary part is an oscillation (similar to the temporal one).

Both real and imaginary parts of this function solve the equation (as any linear combination does). So we take the real part and insert β – only the solution for β with negative sign makes sense as the other one would describe temperature increasing to infinity.

$T(t,x) = Re \left(T_{0}e^{i\omega t} e^{-\sqrt{\frac{\omega}{2D}}(1+i)x}\right)$

The thing in the exponent has to be dimension-less, so we can express the combinations of constants as characteristic lengths, and insert the definition of ω=2π/τ):

$T(t,x) = T_{0} e^{-\frac{x}{l}}cos\left(2\pi\left(\frac {t} {\tau} -\frac{x}{\lambda }\right)\right)$

The two lengths are:

• the wavelength of the oscillation $\lambda = \sqrt{4\pi D\tau }$
• and the attenuation length  $l = \frac{\lambda}{2\pi} = \sqrt{\frac{D\tau}{\pi}}$

So the ratio between those lengths does not depend on the properties of the material and the wavelength is always much shorter than the attenuation length. That’s why there is hardly one period visible in the plots.

The plots have been created with this parameters:

• Heat conductivity κ = 0,0019 kW/mK
• Density ρ = 2000 kg/m3
• Specific heat cp = 1,3 kJ/kgK
• tau = 1 year = 8760 hours

Thus:

• Diffusivity D = 0,002631 m2/h
• Wavelength λ = 17 m
• Attenuation length l = 2,7 m

The wave (any wave) propagates with a speed v equivalent to wavelength over period: v = λ / tau.

$v = \frac{\lambda}{\tau} = \frac{\sqrt{4\pi D\tau}}{\tau} = \sqrt{\frac{4\pi D}{\tau}}$

The speed depends only on the period and the diffusivity.

The maximum of the temperature as observed in a certain depth x is delayed by a time equal x over v. Cross-checking our measurements of the temperature T(30cm) and T(100cm), I would thus expect a delay by 0,7m / (17m/8760h) = 360 h = 15 days which is approximately in agreement with experiments (taking orders of magnitude). Note one thing though: Only the square root of D is needed in calculations, so any error I make in assumptions for D will be generously reduced.

I have not yet included the geothermal linear temperature gradient in the calculation. Again we are grateful for linearity: A linear – zero-curvature – temperature profile that does not change with time is also a trivial solution of the equation that can be added to our special exponential solution.

So the full solution shown in the plot is the sum of:

• The damped oscillation (oscillating about 0°C)
• Plus a constant denoting the true yearly average temperature
• Plus a linear decrease with depth, the linear correction being 0 at the surface to meet the boundary condition.

If there would be no geothermal gradient (thus no flow from beneath) the temperature at infinite distance (practically in 20m) would be the same as the average temperature of the surface.

Daily changes could be taken into account by adding yet another solution that satisfies an amendment to the boundary condition: Daily fluctuations of temperatures would be superimposed to the yearly oscillations. The derivation would be exactly the same, just the period is different by a factor of 365. Since the characteristic lengths go with the square root of the period, yearly and daily lengths differ only by a factor of about 19.

________________________________________

Intro to geothermal energy:

Geothermal gradient and energy of the earth:

These data for bore holes using one scale show the gradient plus the disturbed surface region, with not much of a neutral zone in between.

Theory of Heat Conduction

Heat Transfer Equation on Wikipedia
Textbook on Heat Conduction, available on archive.org in different formats.

I have followed the derivation of temperature waves given in my favorite German physics book on Thermodynamics and Statistics, by my late theoretical physics professor Wilhelm Macke. This page quotes the classic on heat conduction, by Carlslaw and Jäger, plus the main results for the characteristic lengths.

# Ice Storage Challenge: High Score!

Released from ice are brook and river
By the quickening glance of the gracious Spring;
The colors of hope to the valley cling,
And weak old Winter himself must shiver,
Withdrawn to the mountains, a crownless king.

These are the first lines of the English version of a famous German poem on spring, from the drama Faust, by Johann Wolfgang von Goethe. Weird factoid about me: I was once inclined to study literature, rather than physics. But finally physics won, so this is a post about joyful toying with modeling heat transport in ice and water.

After 46 days we had a high score: The ice cube, generated by our heat pump, stopped growing at about 15m3. About 10mof water remained unfrozen. After the volume of ice had been in a steady state for a a while, we turned on the solar collector again to return to standard operations.

Where did the energy for the heat pump come from before?

The lid of the tank is insulated against ambient air, the solar collector was not operational, and no ice had been created: The remaining energy has to be provided by the 5th element that cannot be shut off: 1) water 2) ice, 3) ambient air, 4) solar radiation … 5) ground.

Normally ground supplies about 15 W per m2 surface area – deduced from monitoring the power transported with the brine flow and energy accounting for the tank. The active interface between tank and ground below frost depth is about 35 m2. This results in about 0,5 kW in total, thus just 12 kWh per day, much lower than the ~ 50 kWh ambient energy fed into the heat pump.

After much deliberation and playing with the heat transfer equation we came up with this description of the evolution of the ice cube:

Phase 1: Growth of ice into water.

• Ice starts to grow from the heat exchanger tubes into the remaining water. These tubes are installed in a meandering pattern, traversing the storage tank.
• At some point the thick layers of ice covering adjacent parts of the pipes touch each other. The surface of this solid ice cube is smaller than the interface between the meandering ice formations and water before. The power needed by the heat pump has to be pushed through a smaller surface – which is only possible if the temperature gradient within the ice gets larger. As the temperature at the ice-water interface has to be 0°C, the temperature at the heat exchanger has to decrease. This is exactly what we see from monitoring data – brine temperature drops well below 0°C.
• Side-effect: Due to the lower brine temperature the coefficient of performance  decreases slightly. So more of the total heating energy needs to be provided by the electrical input. We call this the heat source paradox: The worse performance is, the more you spare the energy stored in the heat source. Thanks to this self-protection mechanism, the energy in the tank will not suddenly drop to zero.

Evolution of the volume of ice, ambient temperature and brine temperature over time. The ice is now quickly melting again – in March the collector is already harvesting enough energy again for balancing heating demands!

Phase 2: Ice touching ground.

• As long as there is some water between ice and ground, the water temperature is 0°C. This is the temperature ground ‘sees’ and the temperature which is relevant for the low heat transport from ground to water.
• Ice touches some surfaces of the cuboid tank – the ones where the heat exchanger tubes are closest to the surface. Now ground is directly connected to ice with its temperatures < 0°C. The temperature gradient between ground and ice provides for a higher flow of energy. This is also indicated by the evolution of the temperature in the ground below the tank: While temperatures of undisturbed ground and the region below the tank had been aligned before, ground temperature beneath the tank still kept getting lower – although a few meters away from the tank ground is already warming up again.
• If enough heat is delivered by ground, no more heat is needed by freezing the remaining water in the tank. When ground temperature reaches zero, it can even freeze – which happens with geothermal systems, too. We might have extended the ice storage into ground.

The temperature sensor closest to the tank in 30 cm underneath. A few meters from the tank ground is already warming up again (following the standard ‘yearly temperature wave’), but below the tank the temperature is still getting lower – as highlighted by the blue rectangle.

Heat transport within ice is actually more efficient than transport in water: Ice has 4 times the heat conductivity of water, and 10 times the thermal diffusivity. The latter is a measure for the time a deposited ‘lump of heat’ will be spread in space:

So we have built a very efficient cold bridge between the heat exchanger and ground. Everything is consistent with the poetry of the differential equation of heat transfer.

I marvel at the intriguing and mathematically appealing physics in my backyard!

For the backyard (‘Office desk farming’). Happy Easter everybody!