We report a modeling and numerical simulation study of density-driven natural convection during geological CO2 storage in heterogeneous formations. We consider an aquifer or depleted oilfield overlain by gaseous CO2, where the water density increases due to CO2 dissolution. The heterogeneity of the aquifer is represented by spatial variations of the permeability, generated using Sequential Gaussian Simulation method. The convective motion of the liquid with dissolved CO2 is investigated. Special attention is paid to instability characteristics of the CO2 concentration profiles, variation of mixing length, and average CO2 mass flux as a function of the heterogeneity characterized by the standard deviation and the correlation length of the log-normal permeability fields. The CO2 concentration profiles show different flow patterns of convective mixing such as gravity fingering, channeling, and dispersive based on the heterogeneity medium of the aquifer. The variation of mixing length with dimensionless time shows three separate regimes such as diffusion, convection, and second diffusion. The average CO2 mass flux at the top boundary decreases quickly at early times then it increases, reaching a constant value at later times for various heterogeneity parameters.