Accurate phase behavior modeling of hydrocarbon and aqueous mixtures plays a critical role in simulation of compositional flow in subsurface reservoirs, such as miscible gas flooding and CO2 sequestration. As Michelsen proposed his groundbreaking works in stability test and phase split calculation, PT flash calculation has been well developed in the past four decades and become the most popular flash technique. However, as research interests move to more complicated reservoir fluids, some inherent drawbacks of PT flash formulations show up and recent researches focus on a promising alternative called VT flash calculation. In this thesis, VT flash calculation is used in place of PT flash to model phase behaviors of hydrocarbon and aqueous mixtures. A dynamical model, together with a thermodynamically stable numerical algorithm, is developed to calculate equilibrium phase amounts and compositions with/without capillary effect to simulate phase behaviors of unconventional/conventional hydrocarbon mixtures. In order to model water-containing mixtures, the cubic equation of state is replaced by the Cubic-PlusAssociation equation of state, and a salt-based Cubic-Plus-Association model is developed to calculate phase behaviors of CO2-brine systems. The combination of VT flash calculation and the salt-based Cubic-Plus-Association model accurately estimate CO2 solubility in both single- and mixed-salt solutions, and it exhibits close prediction accuracy with a more sophisticated electrolyte Cubic-Plus-Association model. At the end, the ultimate goal is to develop an efficient two-phase VT-flash compositional flow algorithm. The multilayer nonlinear elimination method is used to remove locally high nonlinearities based on the feedback of intermediate Newton solutions. To further improve the computational efficiency, a modified shadow region method is used to bypass unnecessary stability tests. Although nonlinear elimination fails to fully resolve the convergence issue, which roots in the nondifferentiable equilibrium pressure at the points of phase boundary, the number of time refinements is significantly reduced and the improved VT-flash compositional flow algorithm with multilayer nonlinear elimination method successfully simulates a number of numerical examples with and without gravity.