This paper presents a numerical implementation of the coupled hydro-thermal-mechanical logic in the commercial code FLAC3D. The numerical model uses the Boussinesq approximation in which fluid density variations are neglected in all but the body force term of the equation of motion. The numerical solution is compared to an analytical solution (obtained using stability analysis) for a convection cell in a confined layer heated from below. Numerical simulation examples are presented to illustrate the development of various cell configurations in a layered system heated from below. The numerical model can be applied to the analysis of geothermal groundwater convection in sedimentary basins. The study of the convection mechanism is important because it provides a natural heat exchanger that can be accessed without engineered hydraulic fracturing because of high natural permeability.