In the petroleum industry, accurately estimating wellbore heat loss in thermal recovery processes remains a critical problem. One difficulty lies in simulating heat transfer and fluid dynamics within wellbore annuli. A literature survey shows that the state-of-the-art thermal wellbore simulators use empirical correlations to calculate the heat loss through wellbore annuli. As more sophisticated wells have been drilled, there is a growing need for a more detailed wellbore annulus heat transfer model. In this study, a 2D transient mathematical model is proposed for the conjugate natural convection and radiation within wellbore annuli. The governing equations consist of a continuity equation, a vorticity transfer equation, an energy transport equation and a radiative transfer equation (RTE). A finite volume approach with a second-order upwinding scheme is implemented for discretization. Newton-Raphson iterations are deployed for linearization. The algorithm is validated by consistence in simulation results compared with benchmark numerical solutions with the Rayleigh number up to 107. Parameters such as an aspect ratio, a radius ratio and a conduction-to-radiation coefficient are examined. A case study on vacuum insulated tubing heat transfer using Marlin Well A-6 data demonstrates the merits of the developed program by the consistence of simulation results compared with field measurements.