Vacuum polarization of massive scalar fields in a thermal state at arbitrary temperature is studied near the horizon of a Reissner–Nordström black hole. We derived an analytic form of <ϕ2> approximately in the large mass limit near the black hole horizon. We uses the zeroth order WKB approximation and power series expansion near the horizon for the Euclideanized mode function. Our formula for the vacuum polarization shows regular behavior on the horizon if the temperature of the scalar field is equal to the Hawking temperature of the black hole. The finite part of the vacuum polarization agrees with the result of the DeWitt–Schwinger approximation up to O(m-4) which is the next-to-leading order of the expansion.