This paper describes a new distributed hydrological model, called GEOtop. The model accommodates very complex topography and, besides the water balance, unlike most other hydrological models, integrates all the terms in the surface energy balance equation. GEOtop uses a discretization of the landscape based on digital elevation data. These digital elevation data are preprocessed to allow modeling of the effect of topography on the radiation incident on the surface, both shortwave (including shadowing) and longwave (accounting for the sky view factor). For saturated and unsaturated subsurface flow, GEOtop makes use of a numerical solution of the 3D Richards’ equation in order to properly model, besides the lateral flow, the vertical structure of water content and the suction dynamics. These characteristics are deemed necessary for consistently modeling hillslope processes, initiation of landslides, snowmelt processes, and ecohydrological phenomena as well as discharges during floods and interstorm periods. An accurate treatment of radiation inputs is implemented in order to be able to return surface temperature. The motivation behind the model is to combine the strengths and overcome the weaknesses of flood forecasting and land surface models. The former often include detailed spatial description and lateral fluxes but usually lack appropriate knowledge of the vertical ones. The latter are focused on vertical structure and usually lack spatial structure and prediction of lateral fluxes. Outlines of the processes simulated and the methods used to simulate them are given. A series of applications of the model to the Little Washita basin of Oklahoma using data from the Southern Great Plains 1997 Hydrology Experiment (SGP97) is presented. These show the model’s ability to reproduce the pointwise energy and water balance, showing that just an elementary calibration of a few parameters is needed for an acceptable reproduction of discharge at the outlet, for the prediction of the spatial distribution of soil moisture content, and for the simulation of a full year’s streamflow without additional calibration.

Optimal channel networks (OCN’s) obtained by minimizing the local and global rates of energy expenditure evolve automatically from arbitrary initial conditions to network configurations exhibiting fractal and multifractal statistics indistinguishable from those observed in nature. It is suggested that OCN’s are spatial models of self-organized criticality and that natural fractal structures like river networks may arise as a joint consequence of optimality and randomness.

This paper explores the similarities of digital elevation maps (DEMs) of natural river basins and optimal channel network (OCN) configurations obtained minimizing the total rate of energy expenditure in the system as a whole and in its parts. Striking similarities are observed for natural and optimal networks in their fractal aggregation structure and in certain multifractal structures found to be characteristic of river basins. Our results suggest, upon critical assessment of the reliability of the identification of the attractor of the underlying dynamics implied by our optimality concepts, that fractal structures are indeed possibly a product of least energy dissipation. Power laws emerging in the description of the distribution of aggregated quantities from both DEMs and OCNs suggest a link with the framework of self‐organized criticality in the dynamics of natural channel network formation. Also, the geomorphological description of OCNs reveals surprising analogies with well‐known empirical or experimental results, A comparison of Peano's basins with OCNs suggests that nature seems to reject the type of strict self‐similarity exhibited by Peano's construct in favor of different shapes implying statistical self‐similarity not only because of chance acting through random conditions but also because of necessity as reflected by least energy expenditure considerations.

Phenomena involving frozen soil or rock are important in many natural systems and, as a consequence, there is a great interest in the modeling of their behavior. Few models exist that describe this process for both saturated and unsaturated soil and in conditions of freezing and thawing, as the energy equation shows strongly non-linear characteristics and is often difficult to handle with normal methods of iterative integration. Therefore in this paper we propose a method for solving the energy equation in freezing soil. The solver is linked with the solution of Richards equation, and is able to approximate water movement in unsaturated soils and near the liquid-solid phase transition. A globally-convergent Newton method has been implemented to achieve robust convergence of this scheme. The method is tested by comparison with an analytical solution to the Stefan problem and by comparison with experimental data derived from the literature.