When modelling fluid fl ow in subsurface, the impact of solid deformation on fluid fl ow is often oversimplified/neglected in reservoir simulators. It is assumed that solid volume/s-tate of stress is a function of fluid pressure, while the opposite effect is not considered. This assumption is made mainly to reduce computational costs and complexity of fl uid flow models. Nevertheless, this simplification is not valid in case of unconsolidated rocks. This oversimplification results in wrong estimation in prediction of surface subsidence, earthquakes and fault activation, fluid production and rock permeability values, etc. Numerous reports suggest that neglecting the two-way coupling (i.e. both fluid-to-solid and solid-to- fluid coupling ) has led to disastrous events in many cases. This necessitates modified fl uid models and simulators which take into account the two-way coupled nature of solid deformation and fl uid flow in porous media. Efforts have been made to model this two-way coupled nature properly which can be categorized as follows. There have been attempts to connect commercial softwares to model this problem. They fail due to differences in data structure, different underlying assumptions embedded and due to the increased computational costs. Therefore, it is essential to integrate fluid modelling and solid deformation into a single simulator. This requires developing new numerical models. Presuming an elastic nature for unconsolidated rocks, Biot's consolidation equations are employed to numerically model the two-way coupled solid deformation and fluid fl ow in porous media, so called poroelasticity. During my master work, I developed 2D MATLAB codes based on two different finite volume discretization schemes (cell-centred and vertex-centred FVM). In the first stage, two cell-centred FVM 2D MATLAB codes were developed: One to model fluid flow, and one to model solid equations in poroelasticity. At the next stage, the two cell centred codes were iteratively coupled. Then, another 2D fully coupled model based on a vertex-centred finite volume discretization scheme was developed for poroelasticity. The fl uid and solid data structure in both developed models are the same; in other words, unknowns are collocated. To verify and conduct error analysis on the developed finite volume based simulators, 2D MATLAB codes for one classic benchmark problem in poroelasticity, namely the Mandel problem, as well as 2D MATLAB codes for synthetic test cases were developed. At the next stage, the performance of the two models is compared. Though both methods illustrate adequate performance, fully coupled finite volume method is preferred. Finally, the reaction of the fully coupled finite volume model to different systems, its performance under stress and displacement boundary condition configuration, sensitivity of the model to uncertainty of input parameters, robustness of the model, and applications of this method are being analysed. Furthermore, thorough discussion on enhancement of the model is provided. To our knowledge, it is for the first time that a cell-centred finite volume scheme is applied to poroelastic problems and is compared with the vertex-centred finite volume method. In addition, there is no prior investigation done of robustness for finite volume methods. The results showed that both cell-centred and vertex-centred FVM are computationally efficient in applications for poroelastic problems. However, the fully coupled vertex-centred finite volume model outperformed the sequentially coupled cell-centred model in terms of computational efficiency.