CO2 geo-sequestration (CGS) is considered to be a feasible technology for reducing the amount of CO2 emission into the atmosphere. Selection of an appropriate reservoir is vital and requires appropriate knowledge of the involved phenomena and processes. In a CO2 geo-sequestration process, carbon dioxide goes through mainly four storage (trapping) mechanisms: structural and stratigraphic trapping, residual trapping, solubility trapping and mineral trapping. In this study, focus is placed on modeling the first trapping mechanism, together with corresponding deformation and electrokinetic flow. Multiphase fluid flow due to injection of CO2 in an unsaturated reservoir is accompanied by continuous redistribution of pore pressure and effective stress, causing local and regional deformations and probably major uplifting or subsidence. This flow is also accompanied by electrokinetic flow. In such a system, electrokinetic potentials occur due to the interaction between the formation fluid and the mineral grains. Due to pressure gradients, the flow of the pore fluid produces an advective electric current: such a flow generates an electric field, which produces a counter electric current through the interface, known as the self-potential (SP). Since the electrical conductivity of CO2 is lower than that of the formation brine, it can be detected by measuring the self-potential. Based on this, the SP can be used for monitoring CO2 plume movement, a necessary procedure to ensure that geologic sequestration is both safe and effective. In spite of the versatility of the available numerical tools, attempts to model CO2 geo-sequestration in a region and considering events occurring in local areas lead to enormous demands for computational power. This makes the development of numerical tools for CO2 geo-sequestration not only difficult, but rather expensive. In this study, the governing field equations are derived based on the averaging theory and solved numerically based on a mixed discretization scheme. In this scheme, variables exhibiting different nature are treated using different numerical discretization techniques. Techniques such as the standard Galerkin finite element method (SG), the extended finite element method (XFEM), the level-set method (LS) and the Petrov-Galerkin method (PG) are integrated in a single numerical scheme. SG is utilized to discretize the deformation and the diffusive dominant field equations, and XFEM, together with LS, are utilized to discretize the advective dominant field equations. The level-set method is employed to trace and locate the CO2 plume front, and the XFEM is employed to model the associated high gradient in the saturation field front. The use of XFEM for the advective field leads to a computationally efficient, stable and effectively mesh-independent discretization. However, it gives rise to an extra degree of freedom. The use of SG for the deformation and the diffusive fields requires only standard degrees of freedom, limiting the total number of degrees of freedom and making the scheme computationally efficient. Several verification and numerical examples are presented for both homogenous and fractured reservoirs. The examples demonstrate the capability of the proposed mixed discretization model to simulate challenging, coupled analyses. It has been shown that this model is capable of solving problems, which typically involve several state variables with different transient nature, using relatively coarse meshes.