We employ a three-dimensional indirect boundary element method (BEM) to simulate temperature change around an underground liquefied natural gas storage cavern. The indirect BEM (IBEM) uses fictitious heat source strength on boundary elements as basic variables which are solved from equations of boundary conditions and then used to compute the temperature change at other points in the considered problem domain. The IBEM requires evaluation of singular integration for temperature change due to heat conduction from a constant heat source on a planar (triangular) region. The singularity can be eliminated by a semi-analytical integration scheme. However, it is found that the semi-analytical integration scheme yields sharp temperature gradient for points close to vertices of triangle. This affects the accuracy of heat flux, if they are evaluated by finite difference method at these points. This difficulty can be overcome by a combination of using a direct numerical integration for these points and the semi-analytical scheme for other points distance away from the vertices. The IBEM and the hybrid integration scheme have been verified with an analytic solution and then used to the application of the underground storage.