A new theoretical treatment of phase equilibria is developed and applied to aqueous two-phase systems. It is shown that although osmotic virial expansions, truncated at the second virial coefficient terms, have often been used to rationalise both phase diagrams and protein partitioning behaviours, these are generally unsuitable and higher-order terms are required to represent the many-body interactions which occur in such systems. However, when higher order terms are incorporated in the osmotic virial expansion, it is found that a thermodynamic inconsistency is present and, as a consequence of this, it is proven that the expanded expression cannot be used to solve phase-behaviour problems. As a means of circumventing this difficulty it is suggested that, for practical problems, osmotic virial expansions, using what we have termed 'effective' second virial coefficients, can be used in some situations. The approach proposed for the complete description of the phase behaviour of aqueous two-phase systems consists of the combined application of a recently developed binodal model (Y. Guan, T. H. Lilley and T. E. Treffry, Macromolecules, 1993, 26, 3971) and a second virial coefficient expansion for the solvent. Comparisons of this new approach with experimental phase diagrams have been made and the agreement is satisfactory. Some comments have been made about other approaches to the same problems. The limitations of applying virial expansions to protein distributions are briefly addressed. A statistical mechanical treatment, which links in a simple way, osmotic virial coefficients of the McMillan-Mayer and Hill types, and considers both solvation of the polymers and homotactic interactions of the solvent, is also presented.