Hydrothermal convection in porous media is an essential piece of physics in geothermal reservoirs, and understanding them leads to better development of geothermal energy. We analyze the validity of simulating hydrothermal convection using different formulations of partial differential equations. Using the Elder problem as a benchmark, we found out that the stream function formulation and the velocity formulation are a valid and efficient model of hydrothermal convection. The Nusselt number and entropy production are measurements of the quality of convective heat transfer. The Rayleigh number describes the physical properties of a porous media. We use simulations to investigate further the discrepancy in the Nusselt Rayleigh relationship found in previous experiments. The conclusion is that the multiple steady states of convection pattern in a 3D box are the main reason for the discrepancy found in the Nusselt-Rayleigh relationship.