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 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 use solar energy deposited in summer and stored for winter – I think that ‘geothermal heat pumps’ is a bit of a misnomer.

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.

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/m ^{2}** on continents:

- It is small compared to the energy from the sun: In middle Europe, the sun provides about 1.000 kWh per m
^{2}and year, thus 1.000.000Wh / 8.760h =**144W/m**on average.^{2} - It also much lower than the rule-of-thumb power of ‘flat’ ground loop collectors – about
**20W/m**^{2} - The total ‘cooling power’ of the earth is several 10
^{10}kW: Would the energy not be replenished by radioactive decay, the earth would lose a some seemingly impressive 10^{14}kWh per year, yet this would result only in a**temperature difference of ~10**(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).^{-7}°C

**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)

I start with this:

- 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. - A constant geothermal flow of
**65 mW/m**is superimposed to that.^{2} - 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. 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:

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.

**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.

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. 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:

**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.

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:

The energy density [J/m^{3}] 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/m^{3}] and the specific heat c_{p} [J/kg] of this material. The energy flow per area [W/m^{2}] 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:

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:

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

It oscillates about 0, with an amplitude of T_{0}. 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:

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 β:

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

β 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.

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π/τ):

The two lengths are:

- the wavelength of the oscillation
- and the attenuation length

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/m
^{3} - Specific heat c
_{p}= 1,3 kJ/kgK - tau = 1 year = 8760 hours

Thus:

- Diffusivity D = 0,002631 m
^{2}/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.

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 b**y 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.

________________________________________

Further reading:

**Intro to geothermal energy:**

Where does geothermal energy come from?

**Geothermal gradient and energy of the earth:**

Earth’s heat energy budget

Geothermal gradient

Radius and mass of 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.

Observations, about the first temperature v depth graph

Very interesting that the extreme value point of each graph is crossed by the previous quarter’s graph.

Also, from my more simple minded systems modelling approach, it is clear that these plots are not the steady state plots (impossible to get these !), and that is why they show reisdual oscillations. If as you have assumed (hoped) that your model is correct you can do some plots of damped harmonic functions on top of these (one at a time) and fiddle with the coefficients, which of course will include the starting values, until you get a match. You may get four different sets of matching coefficients, in which case the whole analysis will need a big rethink.

Thanks, Howard! Indeed, the extrema seem to be crossed by the summer/winter curves, I just googled some other random (German) website hat shows the temperature wave, and it really seems to be a feature: http://www.uni-kassel.de/fb14/geohydraulik/Lehre/Num_Mod/Material/Temperaturwelle.png

(From this website of a German on fluid dynamics: http://www.uni-kassel.de/fb14/geohydraulik/Lehre/Num_Mod/NumMod_Sum14.html

There are two coefficients in the solution, the average temperature and the amplitude (or min. and max temperature). I did not try to match those exactly onto the measured curve, as measured temperatures differ every year and the temperature shown here are rather rather high as we just had a record-breaking season – the warmest since meteorological measurements had been recorded. For the plotted theoretical waves, I use temperature values that are slightly lower than than current values and those in the figure based on measurements to show more typical values for our climate, and to compare my curve to the ‘standard textbook curve’ as shown in (German) textbooks about geothermal energy.

Yes, I could change those to match my measured wave, which means just shifting and squeezing the waves a bit, but it would not change the characteristic lengths. The propagation speed of the wave depends only on the thermal diffusivity, not on max/min temperature, so it is justified to compare the speed of the measured wave for a warmer winter to the ‘typical’ wave. (Yes, assuming that diffusivity does not depend that much on the climate in a particular year) – and it is diffusivity I am most interested in.

Of course, in order to confirm this fact (that the simple cosine is really a good model and that D is about constant) I need to vary ambient temperature – which means I need much more and more ‘extreme’ seaons with different temperatures, data I don’t have yet. But, after all, I have used D as an input parameter also for my numerical simulations, and in these simulations I use real-live weather data and thus the boundary temperature is defined. It is a gradual iterative refinement: I am also able to cross-check diffusivity based on my numerical simulations to some extent; I used that D as an input here and tweaked it a bit more.

In the numerical simulation there are lots of other aspects I have had to find pragmatic estimates for, like using an effecitvely 1D model for a 3D situation, so an order of magnitude estimate for D is fine for me. My final criterion is if I can predict the volume of ice in our tank correctly – which depends on heat exchange with ground, but also on many other factors, like all our control logic. So far, the calculations for ice were fine, so I assumed my estimate for D is OK. This ‘model’ for the undisturbed ground served more as an additional confirmation that my D is not totally off-base. But actually, it was just fun to calculate something of some practical relevance analytically only :-)

(Edit : and yes – of course it is a non-steady state solution. I can’t be if the temperature at the boundary isn’t. What could be done is adding the daily, additional damped wave on top of the yearly one – but I am mainly interested in daily averages of temperature.)

When you put a wave signal into a system you get wave signals coming out somewhere or other. How about more data and use the temperature at a specific depth as output. Then a discrete z based model, with some accommodation of the delay could prove interesting.

I think you need a fast ‘signal’ for that, much faster than daily and yearly temperature changes, to that your response can be linked unambiguously to the signal. Something like a thermal response test, done for determining the conductivity of ground (e.g. when building geothermal heat pumps) https://en.wikipedia.org/wiki/Thermal_response_test

But I think, on principle, using the existing natural signal (‘sent’ by the sun / air, with period one year) is the same – if you can wait for one year… and that’s what I did here. Sure, It would be nice to have more data for even larger depths, so that the delay time would be larger, but I will not drill a hole :-)

But it all gets boring(!) after 15 meters, so you don’t need to go much further down. Also, you could get quite a reasonable estimate for the beta from the summer plot, as the earlier part is fairly independent of both the past and thr geothermal thing.

Perhaps a misunderstanding: I have that type of reasonable estimate for beta already (or D or which just differ by factors) – from the last plot, the measured data, by comparing the temporal distance of the maxima to the theoretical value. It is this ‘thermal response test’, for the ‘signal with period 1 year’. You can either measure the damped, spatial distribution, fit the curve to get the characteristic length, or you can measure the temporal development at two points and measure the lag / distance of the maxima. Both give you the same information, a number that depends only on D and on tau=1year.

But it is easier to measure the temporal lag, as you only need two sensors and simply have to wait, whereas you need more sensors to monitor at different depths automatically. I just meant it would be nice to have one more sensor in, say, 5 meters depth, but as I am satisfied with the estimate it is not worth it.

Edit / append: BTW we do have older data for more spatial points, taken manually at different depths for some time. I might analyze those for a future post, some day if I have time. But my intention here was rather to check if those automated measurements at 30cm and 100cm shown above already provide me with an option to cross-check beta / D quickly – and they do. And this calculation is not my only way of cross-checking it, as the ice formation covered by numerical simulations is also impacted by the heat flow from ground.

Manual measurements (moving one Pt1000 sensor manually to different positions in a tube) are time-consuming as you need a minute or two for the temperature at a point to reach steady-state; automated measurements at a new position require a new cable from every sensor somewhere outside to the controller room and a free ‘input slot’ at the control unit.