A nonlinear model is studied which describes the evolution of a landscape under the effects of erosion and regeneration by geologic uplift by mean of a simple differential equation. The equation, already in wide use among geomorphologists and in that context obtained phenomenologically, is here derived by reparametrization invariance arguments and exactly solved in dimension d = 1. Results of numerical simulations in d = 2 show that the model is able to reproduce the critical scaling characterizing landscapes associated with natural river basins. We show that configurations minimizing the rate of energy dissipation (optimal channel networks) are stationary solutions of the equation describing the landscape evolution. Numerical simulations show that a careful annealing of the equation in the presence of additive noise leads to configurations very close to the global minimum of the dissipated energy, characterized by mean field exponents. We further show that if one considers generalized river network configurations in which splitting of the flow (i.e., braiding) and loops are allowed, the minimization of the dissipated energy results in spanning loopless configurations, under the constraints imposed by the continuity equations. This is stated in the form of a general theorem applicable to generic networks, suggesting that other branching structures occurring in nature may possibly arise as optimal structures minimizing a cost function.