A conservative sequential fully implicit method is derived for compositional reservoir simulation. Multi-phase flow in porous media comprises coupled complex processes: i.e. elliptic flow equation, hyperbolic transport equation and highly nonlinear phase equilibrium equation. These processes contain very different mathematical characteristics that cannot be efficiently solved by one numerical method. As a result, the fully implicit method may become numerically complex and inefficient because the Jacobian includes the derivatives w.r.t. the variables from all of the different processes involved. Jenny et al. (2004) [12] and Lee et al. (2015) [20] demonstrated that flow (pressure) and transport (saturation) for multi-phase flow without compositional effect can be efficiently solved by a sequential fully implicit method. However, the characteristics of the phase equilibrium equations are very different from those of the transport equations. This paper proposes an iterative method that solves the flow, transport and phase equilibrium equations in a sequential manner. The transport of hydrocarbons through porous media is governed by the multi-phase Darcy's equation, which is used to compute the phase velocities. The hydrocarbon components belonging to the same phase are transported with the same phase velocity. Upon arrival in the destination grid cell, these components are redistributed via a phase equilibrium calculation. This observation leads to simplification of the governing equations by reducing primary variables to four (i.e., pressure and three phase saturations). The nonlinear solution scheme composed of the stages outlined above is proven to preserve mass conservation, while a new degree of freedom, “thermodynamic flux”, is introduced to ensure volume conservation. The sequential algorithm is solved iteratively until pressure, saturation, and phase composition are fully converged. It is well-known that sequential solution schemes may require many iterations or fail to converge if the phase equilibrium calculation involves phase transition with a large volume change. This indicates that the current governing equations may not adequately describe fluid flux during rapid phase transition. With numerical examples we demonstrate that such numerical difficulties are successfully resolved via the thermodynamic flux term.