Module 0

Mathematical & Physical Foundations

Thermodynamics, transport theory, continuum mechanics, and reaction-diffusion systems as applied to living trees

0.1 Irreversible Thermodynamics

Classical equilibrium thermodynamics describes states of minimal free energy, but living trees are perpetually out of equilibrium—they drive fluxes of water, carbon, nitrogen, and signals across steep gradients. The correct framework is non-equilibrium thermodynamics, pioneered by Lars Onsager and Ilya Prigogine.

Entropy Production

The local rate of entropy production per unit volume, \(\sigma\), measures irreversibility. For a system with multiple simultaneous fluxes:

\[ \sigma = \sum_i J_i X_i \geq 0 \]

where \(J_i\) are generalised fluxes (water flow, heat flow, ion current) and \(X_i\) are the conjugate thermodynamic forces (pressure gradient, temperature gradient, electrochemical potential gradient). The inequality is the second law: entropy can only be produced, never destroyed locally.

Expanding for a tree: the full entropy production includes transpirational water flow\(J_w X_w\), carbon assimilation \(J_C X_C\), ionic transport across membranes\(J_{ion} X_{ion}\), and heat dissipation \(J_q X_q\):

\[ \sigma = J_w \nabla \mu_w + J_C \nabla \mu_C + \sum_{ion} J_{ion} \nabla \tilde{\mu}_{ion} + J_q \nabla \left(\frac{1}{T}\right) \geq 0 \]

Linear Phenomenological Laws

When the system is close to equilibrium (small forces), fluxes are linear functions of all forces:

\[ J_i = \sum_k L_{ik} X_k \]

The matrix \(L_{ik}\) contains the phenomenological coefficients. Diagonal terms \(L_{ii}\)couple each flux to its conjugate force (e.g., Darcy's law, Fourier's law). Off-diagonal terms\(L_{ik}\) describe cross-effects: the Soret effect (temperature gradient drives mass flux), electroosmosis (pressure gradient drives ion current), and thermoosmosis.

Onsager Reciprocal Relations

Onsager (1931) proved from microscopic time-reversal symmetry that the phenomenological matrix is symmetric:

\[ L_{ik} = L_{ki} \]

This is profoundly useful: it reduces the number of independent measurements needed, and constrains cross-coupling in biological membranes. For the water–solute system in a semipermeable membrane:

\[ J_w = L_p \Delta P - L_p \sigma_s \Delta \pi, \quad J_s = C_s(1-\sigma_s) J_w + \omega \Delta \pi \]

where \(L_p\) is the hydraulic conductivity, \(\sigma_s\) is the reflection coefficient (Staverman), \(\omega\) is the solute permeability, and \(\Delta \pi\) is the osmotic pressure difference. The Onsager symmetry gives \(L_{ws} = L_{sw} = -L_p \sigma_s C_s\).

0.2 Diffusion and Transport

Fick's First Law

The molar flux of species \(i\) down its concentration gradient is:

\[ \mathbf{J} = -D \nabla c \]

where \(D\) is the diffusion coefficient (m\(^2\) s\(^{-1}\)). In wood,\(D\) is highly anisotropic: diffusion along the grain (axial direction) is 10–100× faster than radially or tangentially, because cell lumina and pit apertures are oriented axially. The diffusion tensor becomes diagonal in the principal material directions (L, R, T):

\[ \mathbf{D}_{wood} = \begin{pmatrix} D_L & 0 & 0 \\ 0 & D_R & 0 \\ 0 & 0 & D_T \end{pmatrix} \]

Fick's Second Law

Combining Fick's first law with mass conservation \(\partial c / \partial t + \nabla \cdot \mathbf{J} = 0\):

\[ \frac{\partial c}{\partial t} = D \nabla^2 c \]

The fundamental solution (Green's function) in 1D for an instantaneous point source at the origin is:

\[ c(x,t) = \frac{M}{\sqrt{4\pi D t}} \exp\!\left(-\frac{x^2}{4Dt}\right) \]

where \(M\) is the total amount released per unit area. The characteristic diffusion length scales as \(\ell \sim \sqrt{4Dt}\), showing that diffusion is inefficient over long distances—the reason trees invest in bulk flow systems (xylem, phloem).

The Péclet Number

The Péclet number compares advective to diffusive transport:

\[ \mathrm{Pe} = \frac{vL}{D} \]

When \(\mathrm{Pe} \gg 1\), bulk flow dominates and concentration profiles are sharp (convection-dominated). When \(\mathrm{Pe} \ll 1\), diffusion smooths gradients rapidly. In the phloem sieve tube, xylem sap velocity \(v \approx 1\) mm s\(^{-1}\), path length \(L \approx 10\) m,\(D_{sucrose} \approx 5 \times 10^{-10}\) m\(^2\) s\(^{-1}\), giving\(\mathrm{Pe} \approx 2 \times 10^{7}\)—completely advection dominated.

Advection–Diffusion Equation

Including both bulk flow velocity \(\mathbf{v}\) and a source/sink term \(R\):

\[ \frac{\partial c}{\partial t} + \mathbf{v} \cdot \nabla c = D \nabla^2 c + R(c, \mathbf{x}, t) \]

This equation governs carbon transport in phloem, ion transport in xylem sap, and hormone signalling throughout the tree body.

0.3 Continuum Mechanics for Wood

Stress and Strain Tensors

The Cauchy stress tensor \(\boldsymbol{\sigma}\) has components \(\sigma_{ij}\) representing the traction in direction \(j\) acting on a surface with normal in direction \(i\):

\[ \boldsymbol{\sigma} = \begin{pmatrix} \sigma_{LL} & \sigma_{LR} & \sigma_{LT} \\ \sigma_{RL} & \sigma_{RR} & \sigma_{RT} \\ \sigma_{TL} & \sigma_{TR} & \sigma_{TT} \end{pmatrix} \]

in the longitudinal (L), radial (R), and tangential (T) wood directions. Angular momentum balance requires \(\sigma_{ij} = \sigma_{ji}\), so there are 6 independent stress components.

The linearised strain tensor for small deformations from displacement field \(\mathbf{u}\):

\[ \epsilon_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) \]

Orthotropic Elasticity of Wood

Wood is an orthotropic material with three planes of symmetry (LR, LT, RT). The generalised Hooke's law relating the 6-component strain vector to the 6-component stress vector via the compliance matrix \(\mathbf{S}\) (9 independent constants):

\[ \begin{pmatrix} \epsilon_L \\ \epsilon_R \\ \epsilon_T \\ \gamma_{RT} \\ \gamma_{LT} \\ \gamma_{LR} \end{pmatrix} = \begin{pmatrix} \frac{1}{E_L} & -\frac{\nu_{RL}}{E_R} & -\frac{\nu_{TL}}{E_T} & 0 & 0 & 0 \\ -\frac{\nu_{LR}}{E_L} & \frac{1}{E_R} & -\frac{\nu_{TR}}{E_T} & 0 & 0 & 0 \\ -\frac{\nu_{LT}}{E_L} & -\frac{\nu_{RT}}{E_R} & \frac{1}{E_T} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{RT}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{LT}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{LR}} \end{pmatrix} \begin{pmatrix} \sigma_L \\ \sigma_R \\ \sigma_T \\ \tau_{RT} \\ \tau_{LT} \\ \tau_{LR} \end{pmatrix} \]

The 9 independent elastic constants are: 3 Young's moduli (\(E_L, E_R, E_T\)), 3 shear moduli (\(G_{LR}, G_{LT}, G_{RT}\)), and 3 independent Poisson's ratios (\(\nu_{LR}, \nu_{LT}, \nu_{RT}\)). Typical softwood values: \(E_L \approx 12\) GPa, \(E_R \approx 1.5\) GPa,\(E_T \approx 0.8\) GPa. The extraordinary longitudinal stiffness arises from aligned cellulose microfibrils in the dominant S2 cell wall layer.

Stress–Strain Behaviour of Wood

In the elastic region the slope gives Young's modulus \(E_L\). Wood yields when the cell wall buckles (in compression) or when bonds in the middle lamella fail (in tension). Fracture is often catastrophic in tension, but progressive in compression (kinking of fibres).

0.4 Reaction–Diffusion Systems & Turing Instability

Many spatial patterns in tree development—wood ray spacing, stomatal patterning, leaf venation— arise from reaction-diffusion (RD) mechanisms. The canonical two-component system is:

\[ \frac{\partial u}{\partial t} = D_u \nabla^2 u + f(u, v), \qquad \frac{\partial v}{\partial t} = D_v \nabla^2 v + g(u, v) \]

where \(u\) is the activator and \(v\) the inhibitor. Alan Turing (1952) showed that a homogeneous steady state stable to spatially uniform perturbations can be destabilised by diffusion when \(D_v \gg D_u\) (the inhibitor diffuses much faster). This counterintuitive result is the Turing instability.

Linear Stability Analysis

Let \((u_0, v_0)\) be the homogeneous steady state. Perturb:\(u = u_0 + \tilde{u} e^{\lambda t + i\mathbf{q}\cdot\mathbf{x}}\),\(v = v_0 + \tilde{v} e^{\lambda t + i\mathbf{q}\cdot\mathbf{x}}\). The linearised system gives the dispersion relation:

\[ \lambda^2 - \lambda \mathrm{tr}(\mathbf{A}_q) + \det(\mathbf{A}_q) = 0 \]

where \(\mathbf{A}_q = \mathbf{J} - q^2 \mathbf{D}\), \(\mathbf{J} = \begin{pmatrix} f_u & f_v \\ g_u & g_v \end{pmatrix}\) is the Jacobian at steady state, and\(\mathbf{D} = \mathrm{diag}(D_u, D_v)\). A Turing instability requires:

1. \(\mathrm{tr}(\mathbf{J}) = f_u + g_v < 0\) (stable without diffusion)
2. \(\det(\mathbf{J}) = f_u g_v - f_v g_u > 0\) (stable without diffusion)
3. \(D_v f_u + D_u g_v > 0\) (cross-kinetics term dominates)
4. \((D_v f_u + D_u g_v)^2 > 4 D_u D_v \det(\mathbf{J})\) (sufficient disparity in diffusion)

The critical wavenumber at onset is \(q_c^2 = \sqrt{\det(\mathbf{J})/(D_u D_v)}\), giving a spatial wavelength \(\Lambda = 2\pi/q_c\). In cambial development, this mechanism may set the spacing of ray initials, with auxin as the activator and cytokinin as a faster-diffusing inhibitor.

0.5 Python: 1D Diffusion Simulation

We simulate Fick's second law in 1D using an explicit finite-difference scheme. The spatial domain represents the cross-section of a cell wall (0 to 100 nm). A pulse of solute is initialised at the left boundary and diffuses rightward.

Python

script.py73 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.cm as cm

# Parameters
L = 100e-9       # domain length 100 nm
nx = 200         # grid points
dx = L / nx
D = 5e-10        # diffusion coefficient (m^2/s), typical for sucrose
dt = 0.4 * dx**2 / D   # stability: dt < dx^2 / (2D)
t_end = 0.02     # total simulation time (s)
nt = int(t_end / dt)

# Initial condition: Gaussian pulse at left side
x = np.linspace(0, L, nx)
c = np.zeros(nx)
sigma0 = 5e-9
c = np.exp(-((x - 10e-9)**2) / (2 * sigma0**2))

# Time stepping (explicit FTCS scheme)
snapshots = []
times_to_save = np.linspace(0, t_end, 8)
save_idx = 0

for n in range(nt):
    t_current = n * dt
    if save_idx < len(times_to_save) and t_current >= times_to_save[save_idx]:
        snapshots.append((t_current * 1e3, c.copy()))
        save_idx += 1
    # Finite difference update (interior)
    c_new = c.copy()
    c_new[1:-1] = c[1:-1] + D * dt / dx**2 * (c[2:] - 2*c[1:-1] + c[:-2])
    # Neumann BCs (no flux at boundaries)
    c_new[0] = c_new[1]
    c_new[-1] = c_new[-2]
    c = c_new

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
fig.patch.set_facecolor('#0a0f0a')
for ax in (ax1, ax2):
    ax.set_facecolor('#0a0f0a')
    for sp in ax.spines.values():
        sp.set_color('#34d399')
    ax.tick_params(colors='#a7f3d0')

colors = cm.YlGn(np.linspace(0.3, 1.0, len(snapshots)))
for i, (t_ms, snap) in enumerate(snapshots):
    ax1.plot(x * 1e9, snap, color=colors[i], lw=1.8, label=f"{t_ms:.3f} ms")

ax1.set_xlabel("Position (nm)", color="#a7f3d0")
ax1.set_ylabel("Normalised concentration", color="#a7f3d0")
ax1.set_title("1D Diffusion: Concentration Profiles", color="#34d399", fontsize=13)
ax1.legend(fontsize=8, facecolor="#0a0f0a", labelcolor="#a7f3d0", edgecolor="#34d399")

# Mean-square displacement over time
times = np.linspace(dt, t_end, 300)
msd = 2 * D * times
ax2.plot(times * 1e3, msd * 1e18, color="#34d399", lw=2.5)
ax2.set_xlabel("Time (ms)", color="#a7f3d0")
ax2.set_ylabel("Mean-square displacement (nm²)", color="#a7f3d0")
ax2.set_title("MSD = 2Dt (Einstein relation)", color="#34d399", fontsize=13)

plt.tight_layout()
plt.savefig("output.png", dpi=120, facecolor="#0a0f0a")
print("Diffusion simulation complete.")
print(f"  Domain: 0 to {L*1e9:.0f} nm, {nx} grid points, dx = {dx*1e9:.2f} nm")
print(f"  D = {D:.2e} m^2/s, dt = {dt*1e6:.3f} us, total time = {t_end*1e3:.1f} ms")
print(f"  Stability criterion (r = D*dt/dx^2): {D*dt/dx**2:.3f} <= 0.5 (OK)")
diffusion_length = np.sqrt(2 * D * t_end) * 1e9
print(f"  Final diffusion length sqrt(2Dt): {diffusion_length:.1f} nm")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Course Overview Next: Xylem Hydraulics →