Analogous to our canopy example above, in calculating the rate of
heat storage in the ground it is useful to mentally divide the soil
profile into a number of vertical layers and calculate the heat flow and
storage within each layer. What matters then is both the thermal
conductivity of the soil and its heat capacity.
The thermal conductivity of the soil depends on its physical
properties, namely the quartz, mineral and organic content of the soil.
These are assumed to be static properties and depend on soil type. The
inbuilt table soilparams in the micropoint
package is a look-up table for deriving these parameters for a given
soil type, and is used by some of the functions we showcase later. An
important point, however, is that it is also dependent on the fraction
of soil water content, which is not fixed and must itself be modelled.
We can demonstrate this dependence using the
thermalConductivity function in the code below. Here the
volume fraction of quartz and minerals (Vq and Vm respectively) and the
mass fraction of clay (Mc) typical of clay loam are provided as inputs.
The temperature and pressure inputs partially determine conductivity in
the air space, and temperature input is also used to check whether the
water is in the form of ice.
Simple estimates
The rate of heat storage in the ground is one of the most difficult
components of the energy balance to calculate explicitly, so we first
present a simplified model before introducing a full formulation. Here,
the soil is assumed to be infinitely deep with uniform thermal
properties, and the surface temperature is approximated as a sinusoidal
cycle. Under these assumptions, ground heat flux is also sinusoidal, but
phase-shifted by 3 hours and scaled according to the amplitude of
surface temperature, the soil’s thermal properties, and the damping
depth, which determines how far temperature fluctuations penetrate into
the soil.
In the code below, we use solve_wholecanopy to estimate
the temperature of a sparsely vegetated surface, assuming this
approximates ground surface temperature. We initially set the heat
storage term (G) to zero, then iteratively update it and pass it back to
solve_wholecanopy on each step. The effect of including
G is shown in the plot. The first iteration (with
G = 0) is shown in red. As the model iterates, intermediate
solutions are plotted as semi-transparent lines transitioning from red
to blue, and the final converged solution is shown in blue. Including
G shifts the temperature cycle slightly later in the day
and reduces its amplitude.
weather <- micropoint::climdata
tme <- as.POSIXlt(weather$obs_time, tz = "UTC")
sel <- which(tme$mon + 1 == 6 &
tme$mday == 21) # select hours of 21 June
weather <- weather[sel,]
groundp <- micropoint::groundparams # ground parameters
vegp <- createplant_inputs("C3") # vegetation parameters for C3 grass
vegp$h <- 0.05 # vegetation height (m)
vegp$pai <- 0.02 # plant area index
Ca <- micropoint::Cafromyear(2017) # atmospheric CO2 concentration (ppm)
# -------------------------- Initialise variables ------------------------- #
G <- rep(0, 24) # ground heat flux (W m^-2)
error <- 1e99 # difference from previous iteration
Told <- 0 # temperature in previous iteration
itr <- 1 # iteration counter
cols <- adjustcolor(colorRampPalette(c("red", "blue"))(12),
alpha.f = 0.25)
tt <- 0:23 # hours of day
om <- 2 * pi / 24 # angular frequency
while (error > 1e-3) {
# ------------------------ Cycle through each hour --------------------- #
Ts <- 0 # initialise canopy temperature
for (hr in 1:24) {
# ---------------------- Run canopy solver -------------------------- #
can <- solve_wholecanopy(weather, vegp, groundp,
hr,
lat = 49.96807,
long = -5.215668,
zref = 2,
G[hr],
theta = 0.3,
soilrh = 0.85,
Ca)
Ts[hr] <- can$tcanopy # canopy temperature
}
# Plot diurnal temperature on each iteration
if (itr == 1) {
plot(Ts ~ tt, type = "l", col = "red",
ylim = c(5, 55),
xlab = "", ylab = "",
lwd = 2)
} else {
plot(Ts ~ tt, type = "l", col = cols[itr],
ylim = c(5, 55),
xlab = "", ylab = "",
lwd = 2)
}
par(new = TRUE)
# ----------- Calculate decimal hour of peak temperature -------------- #
m <- lm(Ts ~ sin(om * tt) + cos(om * tt))
hrmx <- as.numeric(((pi / 2 -
atan2(coef(m)[3], coef(m)[2])) / om) %% 24)
# ------------------------ Calculate ground flux ----------------------- #
A0 <- (max(Ts) - min(Ts)) / 2 # amplitude of temperature cycle
k <- thermalConductivity(Vq = 0.06, # volumetric quartz fraction
Vm = 0.509, # volumetric mineral fraction
Vo = 0.02, # volumetric organic fraction
Vw = 0.3, # volumetric water fraction
Mc = 0.5422, # clay mass fraction
Tc = Ts, # temperature (°C)
pk = weather$pres)
Ch <- heatCapacity(Vq = 0.06, # volumetric quartz fraction
Vm = 0.509, # volumetric mineral fraction
Vo = 0.02, # volumetric organic fraction
Vw = 0.3, # volumetric water fraction
Tc = 15, # temperature (°C)
pk = 101.3)
ka <- k / Ch # thermal diffusivity
DampDepth <- sqrt(2 * ka * 3600 / om) # damping depth (m)
# Compute ground heat flux
G <- sqrt(2) * A0 *
sin(om * (tt - hrmx) + 3 * pi / 4) / DampDepth
error <- max(abs(Ts - Told)) # convergence error
Told <- Ts
itr <- itr + 1
if (itr > 12) itr <- 12 # ensure colours plot correctly
}
# ---------------------- Plot final converged temperature ----------------- #
plot(Ts ~ tt, type = "l", col = "blue",
ylim = c(5, 55),
xlab = "Hour",
ylab = expression("Temperature (" * degree * "C)"),
lwd = 2)
The damping depth also provides a route to computing sub-surface
temperatures. In the code below, we derive these for depths of 5 cm and
20 cm respectively
# Assumes code above already run
Tav <- mean(Ts) # average temperature (°C)
A0 <- (max(Ts) - min(Ts)) / 2 # amplitude of cycle
# ---------------- Temperatures at 5 cm depth ---------------------------- #
z <- -0.05 # depth below surface (m)
Tz5 <- Tav + A0 * exp(z / DampDepth) *
sin(om * (tt - hrmx) + pi / 2 + z / DampDepth)
# ---------------- Temperatures at 20 cm depth --------------------------- #
z <- -0.2 # depth below surface (m)
Tz20 <- Tav + A0 * exp(z / DampDepth) *
sin(om * (tt - hrmx) + pi / 2 + z / DampDepth)
# ----------------------------- Plot results ----------------------------- #
# Actual surface temperature
plot(Ts ~ tt, type = "l", col = "red",
ylim = c(5, 55),
xlab = "Hour",
ylab = expression("Temperature (" * degree * "C)"),
lwd = 2)
par(new = TRUE)
# Assumed sinusoidal cycle
Tz0 <- Tav + A0 * sin(om * (tt - hrmx) + pi / 2)
plot(Tz0 ~ tt, type = "l", col = "red",
ylim = c(5, 55),
xlab = "", ylab = "",
lwd = 2, lty = 2)
par(new = TRUE)
# 5 cm depth
plot(Tz5 ~ tt, type = "l",
ylim = c(5, 55),
xlab = "", ylab = "",
lwd = 2)
par(new = TRUE)
# 20 cm depth
plot(Tz20 ~ tt, type = "l", col = "blue",
ylim = c(5, 55),
xlab = "", ylab = "",
lwd = 2)
legend("topleft",
legend = c("Surface temperature",
"Sinusoidal approximation",
"5 cm depth",
"20 cm depth"),
col = c("red", "red", "black", "blue"),
lty = c(1, 2, 1, 1),
lwd = 2,
bty = "n")
While the procedure above provides a reasonable way to estimate
ground heat storage and subsurface temperatures, it relies on
assumptions that are clearly limiting. In particular, assuming a
sinusoidal surface temperature cycle is unrealistic, and the assumption
of vertically uniform thermal properties is problematic given their
dependence on soil water content. We therefore turn to a more
comprehensive model.
Soil temperature model
Here we divide the soil into multiple layers. The surface layer has a
non-zero energy balance determined by the absorbed radiation and the
exchange of longwave radiation, sensible and latent heat with the air
above. When this balance is positive the surface warms, and when it is
negative it cools. Heat is then conducted between adjacent soil layers
at a rate determined by the thermal conductivity and the temperature
gradient between them, and the resulting temperature change in each
layer depends on its heat capacity.
At each timestep the temperature profile is updated by solving the
discretised heat equation across all layers simultaneously using a
tridiagonal matrix formulation. This yields a midpoint (partly implicit)
solution that is numerically stable over typical timestep lengths. The
calculation is embedded within an iteration because the surface heat
flux depends on the surface temperature, which is itself part of the
solution. A fixed temperature is imposed at the lower boundary, and the
solution is constrained to avoid non-physical oscillations in the
vertical temperature profile.
To set up the model, we first use the function
createsoilplist. The user specifies the soil type, the
number of layers, and the total depth of the soil profile, extending to
a lower boundary where temperatures are assumed to be stable. The
function then constructs a geometrically increasing sequence of layer
thicknesses, so that deeper layers are progressively thicker, improving
numerical stability. For each layer, typical values are assigned for
bulk density, heat capacity, and the fractions of quartz, mineral, clay,
and organic matter. The user can optionally apply a scaling factor to
control the vertical distribution of organic matter near the surface. By
default, bulk density is assumed to increase slightly with depth. Once
created, any of the default parameters can be modified byu the user.
Available soil types are Sand, Loamy sand,
Sandy loam, Loam, Silt loam,
Sandy clay loam, Clay loam,
Silty clay loam, Sandy clay,
Silty clay and Clay. In the example below the
default parameters for clay loam and leave
surface_organicmu at its default value of 3, meaning the
surface layer is given three times as much soil organic material as the
average for the profile.
# -------------- Generate list of soil model parameters --------------- #
soilp <- createsoilplist(soiltype = "Clay loam",
nlayers = 15, # number of soil layers
totalDepth = 2, # total soil depth (m)
slope = 0, # slope (degrees)
aspect = 180, # aspect (degrees)
surface_organicmu = 3, # surface organic matter depth
FreeDrain = TRUE) # free drainage lower boundary
The final input, FreeDrain, is used by the soil water
model described later. Once the soil parameters have been generated, the
model is initialised using InitSoilheatmod as follows:
# Initial soil temperature sequence: surface to lower boundary layer
Te <- seq(18, 11, length.out = 16)
theta <- rep(0.3, 16) # volumetric soil water fraction
soilheat <- InitSoilheatmod(soilp, Te, theta)
This function creates a list containing the state variables required
to run the model: initial soil temperature and water content for each
layer, plus a fixed temperature at the lower boundary. A sensible choice
for this boundary condition is the mean annual temperature. Each call to
the soil heat model advances the system by one timestep and returns an
updated object of the same form, with revised soil temperatures. The
model also returns an estimate of the rate of ground heat storage.
We are now ready to run the soil heat model. In addition to the soil
parameters and state variables, the model requires the radiation
absorbed at the surface, air temperature, relative humidity, and
atmospheric pressure at a reference height, along with the resistance to
heat transfer between the soil surface and that reference height
(zref).
In the example below, we run the model for short grass over 10 days
in June. The function solve_wholecanopy is used to estimate
the radiation absorbed by the soil surface, while
groundresistance provides the resistance to heat transfer.
The soil heat model (SoilHeatModel) is then applied at each
hourly timestep.
We store the surface temperature and compute the rate of ground heat
storage. Because solve_wholecanopy requires an estimate of
ground heat storage as an input, the calculations are embedded within an
iterative loop. The final code block plots the resulting time series,
illustrating the diurnal cycle in surface temperature and ground heat
flux.
# ---------------- Create driving data from inbuilt datasets -------------- #
weather <- micropoint::climdata
groundp <- micropoint::groundparams
vegp <- createplant_inputs("C3")
vegp$h <- 0.1 # vegetation height (m)
vegp$pai <- 0.15 # total plant area index
boundaryT <- mean(weather$temp) # lower boundary temperature (°C)
# ---------------- Select data for first 10 days in June ----------------- #
tme <- as.POSIXlt(weather$obs_time, tz = "UTC")
sel <- which(tme$mon + 1 == 6 &
tme$mday < 11)
weather <- weather[sel,]
# ---------------- Create parameters for soil heat model ----------------- #
soilp <- createsoilplist(soiltype = "Clay loam",
nlayers = 15, # number of soil layers
totalDepth = 2, # total soil depth (m)
slope = 0, # slope (degrees)
aspect = 180, # aspect (degrees)
surface_organicmu = 3, # surface organic matter depth
FreeDrain = TRUE) # free drainage lower boundary
Te <- seq(weather$temp[1],
boundaryT,
length.out = 16) # initial temperature profile (°C)
theta <- rep(0.3, 16) # volumetric water content
soilheat <- InitSoilheatmod(soilp, Te, theta) # initialise soil heat model
# -------------------------- Initialise variables ------------------------- #
error <- 1e99 # difference between current and previous G
G <- 0 # ground heat flux (W m^-2)
Tsurface <- 0 # ground surface temperature (°C)
Gflux <- rep(G, 24) # ground heat flux for each timestep
for (hr in 1:240) {
while (error > 1e-4) {
# ---------------------- Run canopy model ---------------------------- #
can <- solve_wholecanopy(weather, vegp, groundp,
hr,
49.96807,
-5.215668,
zref = 2,
G,
0.3,
0.8)
# -------- Compute resistance from ground surface to zref ----------- #
rHg <- groundresistance(vegp$h,
vegp$pai,
can$uf,
can$LL,
can$psih,
zref = 2)
# ------------------------- Run soil heat model --------------------- #
soilheat <- SoilHeatModel(soilheat,
soilp,
Rabs = can$Rabs,
Tref = weather$temp[hr],
relhum = weather$relhum[hr],
atmPressure = weather$pres[hr],
rHg)
error <- abs(G - soilheat$Gflux) # convergence error
G <- soilheat$Gflux
}
Gflux[hr] <- G
Tsurface[hr] <- soilheat$Te[1]
soilheat$oldTe <- soilheat$Te # update model state
error <- 1e99 # reset convergence criterion
}
# ----------------------------- Plot time series -------------------------- #
par(mfrow = c(2, 2))
tmec <- as.POSIXct(tme[sel], tz = "UTC")
# Soil surface temperature
plot(Tsurface ~ tmec, type = "l", lwd = 2, col = "red",
xlab = "Date",
ylab = expression("Temperature (" * degree * "C)"))
# Ground heat flux
plot(Gflux ~ tmec, type = "l", lwd = 2, col = "blue",
xlab = "Date",
ylab = expression("Ground heat flux (W " * m^-2 * ")"))
# ---------------- Plot normalised diurnal cycle ------------------------- #
# Values are scaled between each day's minimum and maximum
# to emphasise differences in diurnal shape among days
plotdailycycle <- function(v, clr) {
# Convert hourly vector into day × hour matrix
vd <- matrix(v, ncol = 24, byrow = TRUE)
# Daily minima and maxima
vmn <- rep(apply(vd, 1, min), each = 24)
vmx <- rep(apply(vd, 1, max), each = 24)
# Scale values between 0 and 1
vfrac <- (v - vmn) / (vmx - vmn)
# Rebuild day × hour matrix
vdfrac <- matrix(vfrac, ncol = 24, byrow = TRUE)
for (i in 1:dim(vdfrac)[1]) {
# Plot each daily cycle
v24 <- vdfrac[i,]
if (i < dim(vdfrac)[1]) {
plot(v24, type = "l",
ylim = c(0, 1),
col = clr,
xlab = "", ylab = "")
par(new = TRUE)
} else {
plot(v24, type = "l",
ylim = c(0, 1),
col = clr,
xlab = "Hour",
ylab = "Normalised diurnal cycle")
}
}
}
# Plot daily cycles for surface temperature and ground heat flux
plotdailycycle(Tsurface, rgb(1, 0, 0, 0.5))
plotdailycycle(Gflux, rgb(0, 0, 1, 0.5))
Soil water model
Soil hydrological processes are typically represented by dividing the
soil into a set of vertical layers, analogous to the heat model.
Evaporation removes water from the surface, while transpiration extracts
water throughout the profile according to root distribution. At the same
time, water is redistributed within the soil, with fluxes occurring
between layers. However, the gradient that drives these fluxes is not
the gradient in volumetric water content \(\small \theta\), but the gradient in soil
water potential \(\small \psi_w\),
defined as the pressure of water relative to atmospheric pressure (and
therefore negative in unsaturated soils).
A relationship between \(\small
\theta\) and \(\small \psi_w\)
is therefore required. This relationship is highly nonlinear and
soil-specific, and is typically described using an empirical function.
Two widely used formulations are the Campbell and van Genuchten models.
The Campbell model represents the relationship as a simple power law,
defined by an air-entry potential (\(\small
\psi_e\)) and a shape parameter \(\small b\); air-entry potential is the
potential at which air first begins to enter the largest soil pores as
the soil drains. Below this point, water content declines approximately
as a power function of \(\small
\psi_w\). The van Genuchten model uses a more flexible functional
form, with additional parameters that allow it to represent the observed
relationships across a wider range of conditions. Both approaches define
a continuous mapping between one measure and the other, allowing models
formulated in terms of \(\small
\psi_w\) to be solved while retaining \(\small \theta\) for mass balance and
interpretation.
In the code below, we compare the Campbell and van Genuchten
formulations for a clay loam soil. Each model is used to compute the
corresponding water potential across a range of water contents, and the
resulting retention curves are plotted. This illustrates how the two
approaches differ in the shape of the curve linking moisture status to
pressure, particularly across the drying range.
# ------------ Comparison of Campbell and Van Genuchten --------------------- #
Campbell_params <- Campbellparams("Clay loam")
VG_params <- VGparams("Clay loam")
theta <- seq(VG_params$Smin + 1e-3,
VG_params$Smax,
length.out = 100)
psiw1 <- with(Campbell_params,
PsiFromtheta(theta,
psie,
b,
Smax,
MPa = FALSE))
psiw2 <- with(VG_params,
PsiFromthetaVG(theta,
Smin,
Smax,
psie,
alpha,
n))
# ----------------------------- Plot results ------------------------------ #
par(mar = c(5, 5, 2, 2),
mfrow = c(1, 1))
plot(theta, log10(-psiw1),
type = "l",
lwd = 2,
ylim = c(log10(1e9), log10(1.2)),
yaxt = "n",
xlab = expression(theta ~ "(m"^3 * " m"^-3 * ")"),
ylab = expression(psi[w] ~ "(Pa)"),
cex.lab = 2,
cex.axis = 2)
lines(theta, log10(-psiw2),
col = "red",
lwd = 2)
exps <- 0:9
axis(2,
at = log10(10^exps),
cex.axis = 2,
labels = paste0("-", 10^exps))
legend("bottomright",
legend = c("Campbell", "van Genuchten"),
col = c("black", "red"),
lwd = 2,
bty = "n",
cex = 1.5)
You can see that the two approaches match closely in the mid-range of
theta, but differ when the soil is very wet or very dry.
In the code below, we illustrate how the relationship differs between
contrasting soil types by comparing clay and sand using the van
Genuchten formulation. Water potential is computed across a common range
of water contents for each soil, and the resulting curves are plotted.
The comparison shows that, at the same water content, sand has a much
higher (less negative) potential than clay, meaning water is much less
strongly retained in sand than in clay, owing to its larger pore
sizes.
VG_params1 <- VGparams("Clay")
VG_params2 <- VGparams("Sand")
theta <- seq(0.07, 0.37, length.out = 100) # volumetric water content
psiw1 <- with(VG_params1,
PsiFromthetaVG(theta,
Smin,
Smax,
psie,
alpha,
n))
psiw2 <- with(VG_params2,
PsiFromthetaVG(theta,
Smin,
Smax,
psie,
alpha,
n))
# ----------------------------- Plot results ------------------------------ #
par(mar = c(5, 5, 2, 2))
plot(theta, log10(-psiw1),
type = "l",
lwd = 2,
ylim = c(log10(1e12), log10(1.2)),
yaxt = "n",
xlab = expression(theta ~ "(m"^3 * " m"^-3 * ")"),
ylab = expression(psi[w] ~ "(Pa)"),
cex.lab = 2,
cex.axis = 2)
lines(theta, log10(-psiw2),
col = "red",
lwd = 2)
exps <- 0:12
axis(2,
at = log10(10^exps),
cex.axis = 2,
labels = paste0("-", 10^exps))
legend("bottomright",
legend = c("Clay", "Sand"),
col = c("black", "red"),
lwd = 2,
bty = "n",
cex = 1.5)
Hydraulic conductivity determines how readily water can pass through
soil. Unlike saturated conductivity, which is a fixed property for a
given soil, unsaturated hydraulic conductivity declines sharply as the
soil dries. This occurs because water-filled pathways through the pore
network become progressively reduced, so the soil becomes much less able
to transmit water. Hydraulic conductivity is therefore a key control on
the rate of water redistribution through the profile.
In the code below, hydraulic conductivity is calculated across a
range of water contents using the Campbell and van Genuchten
formulations for the same soil. The resulting curves show how
conductivity declines as the soil dries, and allow comparison of the
different functional forms implied by the two models. Note that as in
the example above, a long-linear scale is used to represent hydrualic
conductivity.
Campbell_params <- Campbellparams("Clay loam")
VG_params <- VGparams("Clay loam")
Ks1 <- with(Campbell_params,
hydraulicConductivityFromTheta(theta,
b, # Campbell shape parameter
Smax, # volumetric water content at saturation
Ksat)) # saturated hydraulic conductivity
Ks2 <- with(VG_params,
hydraulicConductivityFromThetaVG(theta,
Smin, # residual water content
Smax, # volumetric water content at saturation
psie, # air-entry water potential
alpha, # van Genuchten scaling parameter
n, # van Genuchten shape parameter
Ksat)) # saturated hydraulic conductivity
# ----------------------------- Plot results ------------------------------ #
par(mar = c(5, 5, 2, 2))
plot(theta, Ks1,
log = "y",
type = "l",
lwd = 2,
ylim = c(1e-20, 1e-4),
xlab = expression(theta ~ "(m"^3 * " m"^-3 * ")"),
ylab = expression("Hydraulic conductivity (m s"^-1 * ")"),
cex.axis = 2,
cex.lab = 2)
lines(theta, Ks2,
col = "red",
lwd = 2)
legend("bottomright",
legend = c("Campbell", "van Genuchten"),
col = c("black", "red"),
lwd = 2,
bty = "n",
cex = 1.5)
With these relationships defined, we can formulate the soil water
model. As with the heat model, the soil is represented as a set of
vertical layers, and fluxes between layers are calculated based on
differences between them. However, unlike heat transfer, the governing
equations for water are strongly nonlinear, primarily because hydraulic
conductivity changes by many orders of magnitude with moisture status.
This means the system cannot be solved directly. Instead, at each time
step the model uses a Newton–Raphson method, which improves an initial
guess by estimating how much each variable must change for the equations
to be satisfied. This requires knowing how the fluxes respond to small
changes in the state variables. For this reason, the Campbell model is
used here, as its simple form makes these sensitivities straightforward
to evaluate, allowing the solver to be implemented efficiently.
As with any vertical flow model, the behaviour at the base of the
soil profile must be specified. This lower boundary condition determines
whether water can leave the system at depth or is effectively
constrained from doing so. In this model, the boundary condition is
controlled by the FreeDrain argument in the soil parameter
list (created using createsoilplist). If
FreeDrain = TRUE, water is allowed to drain freely out of
the bottom of the profile under gravity, representing a deep,
well-drained soil. If FreeDrain = FALSE, the lower boundary
is treated as saturated, equivalent to the presence of a water table at
depth, which limits drainage from the deepest layer.
The soil water model is run at each time step using the function
SoilWaterModel. The function solves the balance between
four processes. Water can enter the soil from rainfall, be lost from the
surface by evaporation, be extracted from different depths by roots, and
move between adjacent soil layers. The first three act as sources or
sinks, while the last redistributes water internally through the
profile.
In the code below we begin by setting up the inputs and initial
conditions for the model. Weather and ground parameters are taken from
the built-in datasets, and simple vegetation characteristics are
defined. The simulation is then restricted to the first 10 days of June.
A soil profile is created with a specified number of layers and depth,
and initial profiles of temperature and water content are assigned.
These are used to initialise the soil heat and soil water model states
and we set FreeDrain = TRUE. A small indexing step
identifies which layer corresponds to approximately 10 cm depth so that
values at this depth can be extracted later.
The model is then advanced in time using an hourly loop. Within each
hour, a second loop iterates until the surface energy balance has
stabilised. This inner loop ensures that the canopy, soil heat, and soil
water components are consistent with one another before moving to the
next time step. Within each iteration, the canopy model is solved first.
This provides canopy temperature, radiation absorption, and the
resistances that control exchanges with the atmosphere. These outputs
are then used to calculate the resistance between the soil surface and
the reference height, which determines how efficiently heat and water
vapour are transferred away from the surface. The soil heat model is
then updated using these conditions. This produces a new temperature
profile and an updated ground heat flux. The difference between
successive estimates of ground heat flux is used to assess convergence,
and the loop continues until this difference becomes small. Once the
energy balance is consistent, the code computes the water fluxes that
drive the soil water model. Transpiration is calculated from the canopy
state using vapour pressure gradients and resistances, and is set to
zero when stomatal resistance is very high. Evaporation is calculated
separately from the top soil layer, using its current water content and
temperature. Together with precipitation from the weather data, these
define the gains and losses of water for that time step.
These fluxes are then passed to the soil water model, along with the
current soil state. The soil temperature profile is also provided so
that temperature-dependent processes are evaluated consistently. The
soil water model is then run, updating water potential and water content
in each layer by balancing evaporation, transpiration, precipitation,
and vertical redistribution. After the soil water model has updated the
profile, the new water contents are passed back to the heat model. This
ensures that the thermal calculations reflect the current water status
of the soil, maintaining the coupling between heat and water.
Once the inner loop has converged, selected outputs are stored for
that hour. In this example, temperature and water content are recorded
at the surface and at approximately 10 cm depth. The current state is
then saved as the previous state, so that the next time step starts from
the updated conditions. After all time steps have been completed, the
stored values are plotted. The results show how temperature and water
content vary through time at the surface and at depth, with surface
conditions responding more rapidly to atmospheric forcing and deeper
layers showing slower, damped responses.
# ---------------- Load example forcing and parameter datasets ------------ #
weather <- micropoint::climdata # meteorological forcing data
groundp <- micropoint::groundparams # ground surface parameters
vegp <- createplant_inputs("C3") # vegetation parameters for C3 grass
vegp$h <- 0.1 # canopy height (m)
vegp$pai <- 0.15 # plant area index
boundaryT <- mean(weather$temp) # lower soil boundary temperature (°C)
# ---------------- Select first 10 days of June -------------------------- #
tme <- as.POSIXlt(weather$obs_time, tz = "UTC")
sel <- which(tme$mon + 1 == 6 &
tme$mday < 11)
weather <- weather[sel,]
# ---------------------- Create soil parameter set ----------------------- #
soilp <- createsoilplist(soiltype = "Clay loam",
nlayers = 15, # number of soil layers
totalDepth = 2, # total soil depth (m)
slope = 0, # slope (degrees)
aspect = 180, # aspect (degrees)
surface_organicmu = 1, # surface organic matter depth
FreeDrain = TRUE) # free drainage lower boundary
n <- length(soilp$b) # number of soil layers
# ---------------------- Initialise soil profiles ------------------------ #
Te <- seq(weather$temp[1],
boundaryT,
length.out = n) # initial temperature profile (°C)
theta <- seq(0.2,
soilp$Smax[n],
length.out = n) # initial soil water profile
soilheat <- InitSoilheatmod(soilp,
Te,
theta) # initialise soil heat model
soilwater <- InitSoilwatermod(soilp,
soilheat,
rootskew = vegp$rootskew) # initialise soil water model
# Identify layer closest to 10 cm depth
nd <- which(soilheat$z[1:10] < 0.1)
nd <- nd[length(nd)] + 1
# ------------------------- Initialise outputs --------------------------- #
Tsurface <- T10cm <- theta_surface <- theta_10cm <- Gflux <- er <- rep(NA, 240)
G <- 0 # initial ground heat flux estimate
error <- 1e99 # convergence criterion
# --------------------------- Main simulation loop ----------------------- #
for (hr in 1:240) {
itr <- 1
# Iteratively solve canopy, soil heat, and soil water balance
while (error > 1e-2 && itr < 20) {
# ---------------------- Solve canopy energy balance ---------------- #
can <- solve_wholecanopy(weather,
vegp,
groundp,
hr,
49.96807,
-5.215668,
zref = 2,
G,
0.3,
0.8)
# -------- Aerodynamic resistance between soil and atmosphere ------- #
rHg <- groundresistance(vegp$h,
vegp$pai,
can$uf,
can$LL,
can$psih,
zref = 2)
# ----------------------- Solve soil heat transport ----------------- #
soilheat <- SoilHeatModel(soilheat,
soilp,
Rabs = can$Rabs, # absorbed radiation (W m^-2)
Tref = weather$temp[hr], # air temperature (°C)
relhum = weather$relhum[hr], # relative humidity (%)
atmPressure = weather$pres[hr], # pressure (kPa)
rHg)
error <- abs(G - soilheat$Gflux) # convergence error
G <- soilheat$Gflux
# ----------------------- Compute canopy transpiration -------------- #
Tk <- can$tcanopy + 273.15 # canopy temperature (K)
if (can$rSl > 1e9) {
# Stomata effectively closed
Transp <- 0
} else {
# Total vapour resistance
rV <- can$rHa + can$rSl
# Saturation vapour pressure inside canopy
es <- satvap(can$tcanopy) * 1000
# Atmospheric vapour pressure
ea <- satvap(weather$temp[hr]) *
(weather$relhum[hr] / 100) * 1000
# Convert vapour flux to mm h^-1
Transp <- (0.018015 / (8.314 * Tk * rV)) *
(es - ea) * 3600
# Prevent negative transpiration
Transp[Transp < 0] <- 0
}
# ----------------------- Compute bare soil evaporation ------------ #
theta_av <- (soilwater$theta[1] +
soilwater$oldtheta[1]) / 2 # average surface water content
Tc_av <- (soilheat$Te[1] +
soilheat$oldTe[1]) / 2 # average surface temperature
Evapmmhr <- evaporation_flux(theta_av,
Tc_av,
weather$temp[hr],
weather$relhum[hr],
rHg,
soilp$psie[1], # air-entry water potential
soilp$b[1], # Campbell shape parameter
soilp$Smax[1], # saturation water content
3600) # timestep (s)
# ----------------------- Assemble forcing data --------------------- #
climdata <- list(Tair = weather$temp[hr], # air temperature (°C)
relhum = weather$relhum[hr], # relative humidity (%)
rHa = rHg, # aerodynamic resistance (s m^-1)
precip = weather$precip[hr], # precipitation (mm h^-1)
Et = Transp, # transpiration (mm h^-1)
Evapmmhr = Evapmmhr) # evaporation (mm h^-1)
# ----------------------- Solve soil water movement ----------------- #
soilwater$Tc <- soilheat$Te
soilwater$oldTc <- soilheat$oldTe
soilwater <- SoilWaterModel(soilwater,
soilp,
climdata,
pTAW = vegp$pTAW, # fraction extractable water
maxNrIterations = 20)
# Update thermal properties using new soil moisture
soilheat$wc <- soilwater$theta
itr <- itr + 1
}
# ---------------------------- Store outputs -------------------------- #
Gflux[hr] <- G
Tsurface[hr] <- soilheat$Te[1]
T10cm[hr] <- soilheat$Te[nd]
theta_surface[hr] <- soilwater$theta[1]
theta_10cm[hr] <- soilwater$theta[nd]
soilheat$oldTe <- soilheat$Te
soilwater$oldvapor <- soilwater$vapor
soilwater$oldtheta <- soilwater$theta
er[hr] <- error
error <- 1e99
}
# --------------------------- Plot time series --------------------------- #
par(mfrow = c(2, 1))
tmec <- as.POSIXct(tme[sel], tz = "UTC")
# Soil temperature time series
plot(Tsurface ~ tmec, type = "l",
lwd = 2,
col = "red",
ylim = c(8, 50),
xlab = "Date",
ylab = expression("Temperature (" * degree * "C)"))
par(new = TRUE)
plot(T10cm ~ tmec, type = "l",
lwd = 2,
col = "blue",
ylim = c(8, 50),
xlab = "",
ylab = "")
legend("topright",
legend = c("Surface", "10 cm"),
col = c("red", "blue"),
lwd = 2,
bty = "n")
# Soil moisture time series
plot(theta_surface ~ tmec, type = "l",
lwd = 2,
col = "red",
ylim = c(0, 0.5),
xlab = "Date",
ylab = expression(theta ~ "(m"^3 * " m"^-3 * ")"))
par(new = TRUE)
plot(theta_10cm ~ tmec, type = "l",
lwd = 2,
col = "blue",
ylim = c(0, 0.5),
xlab = "",
ylab = "")
legend("bottomright",
legend = c("Surface", "10 cm"),
col = c("red", "blue"),
lwd = 2,
bty = "n")
The resulting time series show contrasting behaviour between the
surface and deeper soil layers. Surface temperature responds rapidly to
atmospheric forcing, with large diurnal fluctuations, while temperature
at 10 cm depth is more buffered and varies more slowly. Periods when the
surface soil is wet show noticeably reduced temperature variability,
reflecting the strong buffering effect of higher water content on
thermal properties. A similar pattern is evident for water content. The
surface layer shows sharp declines associated with evaporation, followed
by recovery during wetter periods, whereas the deeper layer remains
comparatively stable. Once the surface becomes wet, it also tends to
remain wet for some time, because hydraulic conductivity is much higher
under these conditions, allowing water to be more readily supplied from
below.
The reader may wish to use the
In the code below we set the air
temperature and pressure at the reference height, assume a sensible heat
flux of 200 W per square metre, and calculate the volumetric heat
capacity of air. We then compute the resistance to zref for each canopy
and use this to estimate its heat-exchange surface temperature. Next, we
obtain resistances for each height vector and, together with the surface
temperatures, derive temperature–height profiles, which are then
plotted.
In the examples below, we show how the
diabatic coefficients are derived iteratively and demonstrate their
influence on the temperature and wind profiles. In contrast to the
previous examples, we set the radiation absorbed the canopy
(
