We show that landscapes and their drainage networks evolve under erosion rules to obey a variational principle. Starting from hydrodynamics we first find by separation of variables a general relationship between erosion and a modified slope-area law. This law encompasses previous work and observations, and we confirm it by simulation. Secondly we show that this corresponds to a solution of a variational minimization model, generalized from Rodriguez-Iturbe et al. Thus a mechanism by which the local process of erosion acts to globally optimize a river network is produced. [S0031-9007 (96)00098-1] PACS numbers: 64.60.Ak, 92.40.Fb Empirical observations on natural rivers by Leopold and Maddock [1] and Leopold [2] revealed power-law relationships between the slope (s), width (w), depth (d), velocity (y), and discharge (Q) of a channel, measured at the same time for different locations amongst river networks in the United States: s~Q 20.5 , w~Q 0.5 , d~Q 0.4 , y~Q 0.1 .The first of these is an observation of the slope-area law. A theoretical framework was put forward by RodriguezIturbe et al. [3] that was able to account for some of these observations. Using two unproven principles an expression was derived for the total rate of expenditure of the mechanical potential energy of all the water in the network (E):The unsubstantiated suggestion was then made that natural river networks evolve in such a way as to make E a global minimum. Computed networks that obey this rule are called optimal channel networks, and these networks have been found to obey Hack's law, Horton's laws, Moon's law, Melton's law, Strahler ordering, and the drainage-basin area law to the same extent as natural networks [4][5][6][7][8]. Different global minimization principles were also suggested by Yang, Song, and Davies and Sutherland [9-11], but in each case a full justification for the method proved difficult to obtain. Rinaldo et al. [12,13], and latterly Sun, Meakin, and Jossang [14], performed a direct minimization of E by computer simulation, using a simulated annealing approach. In this model, a channel network was placed on a square lattice such that the water present at each site in the lattice could drain through the network to the perimeter of the lattice. This network was then rearranged, until a drainage network with a minimum value of E was found. This value was remarkably constant, using a wide variety of initial conditions. The identification s~Q 20.5 was then made [15], enabling a value of height to be allotted to each site in the lattice, and it was discovered that the scalar field of heights was continuous over basin boundaries.In this Letter basic hydrodynamics is used to establish an erosion equation, which is numerically shown to erode the landscape to give steady drainage networks. We then show analytically that this implies a particular generalization of the slope-area law. This law is then revealed to be equivalent to the minimization of a global quantity L, which explains how local erosion can act to minimize a global qu...