A robust solution method for the large-scale, structured sparse nonlinear algebraic equation systems involved in steady-state process modelling calculations is presented. The method retains rapid local convergence by embedding Newton’s method within a branch and bound framework providing rigorous global convergence. The strategy invokes progressively more involved computations as needed, solving linear programs that combine local Jacobian and global bounding information to generate corrected search directions and assess region feasibility. Interval analysis techniques are employed to automatically generate bounding functions from the analytical expressions for the equations. Adaptive domain partitioning and a branch and bound structure provide a formal backtracking mechanism. Details of the algorithm are described, and the results of tests on both small problems and flowsheeting systems in a few thousand variables are presented. The tests demonstrate effective performance on a variety of process modelling problems. In large-scale problems, the importance of using efficient LP solution methods is evidenced.