Vegetation Pattern Formation via Energy-Balance-Constrained Modeling

Abstract

Vegetation in semi-arid environments self-organizes into striking spatial patterns -- bands, spots, labyrinths, and gaps -- with characteristic wavelengths on the order of tens to hundreds of meters. Existing reaction-diffusion models postulate nonlinearities and transport laws from qualitative physical reasoning, making it hard to distinguish essential structural features from artifacts of the chosen forms. Here we show how energy-balance and water-conservation principles can constrain the admissible model class before a specific closure is chosen. These constraints motivate a family of semilinear closures; an Euler--Lagrange representative yields a fourth-order vegetation equation coupled to quasi-steady water transport on a one-dimensional hillslope. Linear stability analysis identifies three instability mechanisms: classical water-mediated feedback, energy-balance spatial coupling, and water deflection by vegetation gradients. Their balance depends on terrain geometry. On slopes, the water-mediated coupling dominates and the model reproduces two empirical observations: pattern wavelength increases with aridity, and vegetation bands migrate uphill. On flat terrain, the energy-balance spatial coupling can drive instability independently. Numerical simulations confirm the linear predictions, and exploratory continuation reveals a narrow hysteresis region consistent with subcritical bifurcation.

Figures (6)

• Figure 1: Dispersion relation $\operatorname{Re}[\sigma(k)]$ (black) and the $k$-dependent increments of each term in the decomposition \ref{['eq:dispersion']}, relative to $k=0$: water-coupling increment $\operatorname{Re}[\mathcal{B}(k)\Phi(k)] - \operatorname{Re}[\mathcal{B}(0)\Phi(0)]$ (blue), energy-balance spatial coupling $Dk^2$ (orange), and fourth-order cutoff $Ek^4$ (pink). Filled circles mark the most unstable wavenumber $k_*$. Both panels use $\Lambda_1 = 0$ and $\Delta = 0$ with other parameters at baseline values from \ref{['tab:baseline']}, including $\nu = 1$. (a) Weak spatial coupling ($\Lambda_2 = 0.005$): $k_* \approx 9.8$, $\sigma_{\max} \approx 0.010$. (b) Strong energy-balance coupling ($\Lambda_2 = 0.1$): the larger interaction length shifts $k_*$ down to $\approx 2.2$ and selects a longer wavelength, consistent with \ref{['eq:kstar']}. These panels use an illustrative regime with $\Lambda_1 = 0$ (suppressing asymmetric energy interactions while retaining water advection), where the local polynomial is the primary determinant of the selected wavenumber. Note the different $k$-axis range in (b).
• Figure 1: Wavelength selection and rainfall dependence on a slope ($\Lambda_1 = 0.10$). (a) Dominant pattern wavelength $\lambda_*$ versus dimensionless rainfall $\rho$: linear theory (blue curve), box-quantized admissible wavelengths $L/n$ in the finite periodic domain (gray markers), and nonlinear PDE simulations (orange points). (b) Dispersion relation $\operatorname{Re}[\sigma(k)]$ at four rainfall values. Filled circles mark the most unstable wavenumber $k_*(\rho)$. As aridity increases (decreasing $\rho$), the most unstable wavenumber $k_*$ shifts to lower values, producing longer wavelengths at lower rainfall consistent with the Deblauwe trend. Note that the $\rho = 0.65$ curve peaks below zero (the uniform state is linearly stable at this rainfall), but the spectral peak still tracks the wavelength that would be selected by subcritical nucleation. Other parameters at baseline values (\ref{['tab:baseline']}).
• Figure 1: Normalized dispersion relations: $\sigma/\sigma_{\max}$ versus $k/k_*$ for the energy balance model (blue solid; $\Lambda_1 = 0$, $\Lambda_2 = 0.01$, other parameters at baseline from Table 5 of the main text including $\nu = 1$) and the Klausmeier model (orange dashed; $a = 1.5$, $m = 0.45$, $v = 182.5$, $D_u = 1$). Both models are evaluated on the upper vegetated branch; the Klausmeier parameters are representative values from the literature Sherratt2005. The energy balance instability band is a sharp spike at $k/k_* = 1$; the Klausmeier band is a broad dome. The normalization by $\sigma_{\max}$ shows only the positive part of each curve.
• Figure 2: Mechanism classification across parameter space. At each point, the full dispersion relation and two reductions---no energy-balance spatial coupling ($\eta = 0$) and no water-mediated feedback ($\beta = \chi = \Delta = 0$)---are evaluated; a mechanism is necessary if removing it eliminates the instability. (a) $(\rho, \Gamma)$ plane. (b) $(\rho, \beta)$ plane. Colors: energy-balance driven (blue), water driven (orange), synergistic (pink), either suffices (light blue), stable (gray). Black contour marks the linear stability boundary. Computed at $\Lambda_1 = 0$ with $\nu = 1$ and other parameters at baseline values (\ref{['tab:baseline']}) unless varied on the axes. This is an illustrative mathematical regime that suppresses the asymmetric energy interaction while retaining water advection; the physical flat-terrain limit is $\Lambda_1 = 0$, $\nu = 0$.
• Figure 2: Spacetime plots of $U(X,T)$ showing band migration. (a) Hillslope ($\nu = 1$, $\Lambda_1 = 0.10$, $\chi = 2$, $\rho = 0.8$): diagonal stripes indicate steady uphill migration. (b) Flat terrain ($\nu = 0$, $\Lambda_1 = 0$, $\chi = 0$, $\Delta = 0$, $\rho = 0.5$): stationary vertical bands. Other parameters at baseline values (\ref{['tab:baseline']}).