This thesis presents mathematical models for two-phase pipeline flow, with an emphasis on applications to CO2 pipeline flow, as well as numerical methods suitable for solving these models. The considered models form a hierarchy of homogeneous (single-velocity) two-phase flow models with relaxation terms that account for transfer processes between the two phases. The relaxation terms model heat, mass and volume transfer caused by differences in temperature, chemical potential and pressure, respectively. The basis of the model hierarchy is a six-equation model with all three relaxation processes present. The rest of the hierarchy is then derived by assuming that one or more of the relaxation processes are infinitely rapid, which results in equilibrium in pressure, temperature and/or chemical potential, which makes a total of eight models. The models are formulated using conservation laws for mass, momentum and energy as well as an advection equation for the gas volume fraction. It is shown that the subcharacteristic condition, which is related to the stability of such models, translates to the requirement that the speed of sound is reduced when a new equilibrium condition is introduced. Expressions for the speeds of sound in the eight models are derived and proven to satisfy the subcharacteristic condition. A mass-transfer model for pipeline flow based on statistical rate theory is derived and formulated as a chemical-potential relaxation term in the pressure-temperature equilibrium model of the hierarchy. The model is used to simulate depressurization of a CO2 pipeline, and the results are found to be quite close to those of the full-equilibrium model. An exponential time-differencing scheme tailor-made for relaxation terms is applied to the model and compared to the Backward Euler method. The exponential time-differencing scheme is an explicit method, but it relies on knowledge of the equilibrium of the relaxation process. The mass-transfer equilibrium value has to be calculated using a Newton-Raphson iteration, which essentially makes both methods implicit, and comparable in both computational cost and accuracy. Finally, the Rankine-Hugoniot-Riemann (RHR) solver is presented, which aims to solve multidimensional conservation laws with source terms. The solver introduces the novel idea of treating flux gradients in other dimensions as additional source terms. The source term and cross-flux term is placed as a singular source in the centre of each cell, which causes a jump in the solution according to a Rankine-Hugoniot condition. The states on either side of a cell interface then define a Riemann problem that is solved by an approximate Riemann solver. The RHR solver is shown to be of second order in space for a 2D scalar advection equation, the 2D isothermal Euler equations and the 2D shallow water equations
