Processes that take place in the Earth’s mantle are not accessible to direct observation. Natural samples of mantle material that have been transported to the surface as xenoliths provide useful information on phase relations and compositions of phases at the pressure and temperature conditions of each rock fragment. In the past, considerable effort has been devoted by petrologists to investigate upper mantle processes experimentally. Results of high temperatures, high pressure experiments have provided insight into lower crust-upper mantle phase relations as a function of temperature, pressure and composition. However, the attainment of equilibrium in these experiments, especially in complex systems, may be very difficult to test rigorously. Furthermore, experimental results may also require extrapolation to different pressures, temperatures or bulk compositions. More recently, thermodynamic modeling has proved to be a very powerful approach to this problem, allowing the deciphering the physicochemical conditions at which mantle processes occur. On the other hand, a comprehensive thermodynamic model to investigate lower crust-upper mantle phase assemblages in complex systems does not exist. ¶ In this study, a new thermodynamic model to describe phase equilibria between silicate and/or oxide crystalline phases has been derived. For every solution phase the molar Gibbs free energy is given by the sum of contributions from the energy of the end-members, ideal mixing on sites, and excess site mixing terms. It is here argued that the end-member term of the Gibbs free energy for complex solid solution phases (e.g. pyroxene, spinel) has not previously been treated in the most appropriate manner. As an example, the correct expression of this term for a pyroxene solution in a general (Na-Ca-Mg-Fe2+-Al-Cr-Fe3+-Si-Ti) system is presented and the principle underlying its formulation for any complex solution phase is elucidated.¶ Based on the thermodynamic model an algorithm to compute lower crust-upper mantle phase equilibria for subsolidus mineral assemblages as a function of composition, temperature and pressure has been developed. Included in the algorithm is a new way to represent the total Gibbs free energy for any multi-phase complex system. At any given temperature and pressure a closed multi-phase system is at its equilibrium condition when the chemical composition of the phases present in the system and the number of moles of each are such that the Gibbs free energy of the system reaches its minimum value. From a mathematical point of view, the determination of equilibrium phase assemblages can, in short, be defined as a constrained minimization problem. To solve the Gibbs free energy minimization problem a ‘Feasible Iterate Sequential Quadratic Programming’ method (FSQP) is employed. The system’s Gibbs free energy is minimized under several different linear and non-linear constraints. The algorithm, coded as a highly flexible FORTRAN computer program (named ‘Gib’), has been set up, at the moment, to perform equilibrium calculations in NaO-CaO-MgO-FeO-Al2O3-Cr2O3-Fe2O3- SiO2-TiO2 systems. However, the program is designed in a way that any other oxide component could be easily added.¶ To accurately forward model phase equilibria compositions using ‘Gib’, a precise estimation of the thermodynamic data for mineral end-members and of the solution parameters that will be adopted in the computation is needed. As a result, the value of these parameters had to be derived/refined for every solution phase in the investigated systems. A computer program (called ‘GibInv’) has been set up, and its implementation is here described in detail, that allows the simultaneous refinement of any of the end-member and mixing parameters. Derivation of internally consistent thermodynamic data is obtained by making use of the Bayesian technique. The program, after being successfully tested in a synthetic case, is initially applied to pyroxene assemblages in the system CaO-MgO-FeO-Al2O3-SiO2 (i.e. CMFAS) and in its constituent subsystems. Preliminary results are presented.¶ The new thermodynamic model is then applied to assemblages of Ca-Mg-Fe olivines and to assemblages of coexisting pyroxenes (orthopyroxene, low Ca- and high Ca clinopyroxene; two or three depending on T-P-bulk composition conditions), in CMFAS system and subsystems. Olivine and pyroxene solid solution and end-member parameters are refined, in part using ‘GibInv’ and in part on a ‘trial and error’ basis, and, when necessary, new parameters are derived. Olivine/pyroxene phase relations within such systems and their subsystems are calculated over a wide range of temperatures and pressures and compare very favorably with experimental constraints.