The task of finding all physically relevant solutions to mathematical models of physical systems remains an important and challenging area of active research in many branches of science and engineering. While there are several useful 'global' methods for finding one or more solutions to models for chemical and other processes, it is the solutions that are generally the primary focus. Singular points are often viewed as something to be avoided and no use is made of the natural connectedness that exists between solutions and singular points. This paper describes a completely different, novel and general approach to finding all physically meaningful solutions and singular points to mathematical models of physical systems by intelligently moving up and down the landscape of the least-squares function. Theoretical foundation for this work rests on the fundamental observations that (1) solutions and singular points are smoothly connected when the model functions are smooth; (2) valleys, ridges, ledges, etc. provide a natural characterization of this connectedness; (3) valleys, ridges, etc. can, in turn, be characterized as a collection of constrained extrema over a set of level curves; (4) there is an equivalent characterization of valleys, ridges, etc. as solutions to generalized, constrained eigenvalue-eigenvector problems; and (5) the natural flow of Newton-like vector fields tends to be along these distinct features of the landscape. Differential geometry is used to provide theoretical support for these fundamental observations and related issues. These observations also form the basis for a new family of algorithms for finding all physically meaningful solutions and singular points called 'global terrain methods', which consist of a series of downhill, equation-solving computations and uphill, predictor-corrector calculations. Downhill movement to either a singular point or solution is conducted using reliable, norm-reducing (complex domain) trust region methods. Uphill movement, on the other hand, is necessarily to a singular point and uses approximate uphill Newton-like predictor steps combined with intermittent corrector steps. Each corrector step is defined by calculating an extremum in the gradient norm on the current level set for the least-squares function, can be shown to be equivalent to a solution to a generalized, constrained eigenvalue-eigenvector problem and helps ensure that valleys and ridges are tracked as closely as desired. Initial starting points are arbitrary, while starting points for subsequent subproblems defining movement from one stationary point to another are along appropriately determined eigendirections, since valleys and ridges are generalized eigenpathways. Collisions with boundaries of the feasible region and the presence of points at infinity are also addressed and the heuristic termination criterion based on the concept of limited connectedness is presented. A variety of numerical results and geometric illustrations for two-dimensional chemical process models are used to make clear key theoretical concepts, to demonstrate the reliability and efficiency of global terrain methods on small scale problems and to show their potential promise on large scale process models problems.