We study the thermal structure of neutron stars with magnetized envelopes composed of accreted material, using updated thermal conductivities of plasmas in quantizing magnetic fields, as well as equation of state and radiative opacities for partially ionized hydrogen in strong magnetic fields. The relation between the internal and local surface temperatures is calculated and fitted by an analytic function of the internal temperature, magnetic field strength, angle between the field lines and the normal to the surface, surface gravity, and the mass of the accreted material. The luminosity of a neutron star with a dipole magnetic field is calculated for various values of the accreted mass, internal temperature, and magnetic field strength. Using these results, we simulate cooling of superfluid neutron stars with magnetized accreted envelopes. We consider slow and fast cooling regimes, paying special attention to very slow cooling of low-mass superfluid neutron stars. In the latter case, the cooling is strongly affected by the combined effect of magnetized accreted envelopes and neutron superfluidity in the stellar crust. Our results are important for interpretation of observations of isolated neutron stars hottest for their age, such as RX J0822-43 and PSR B1055-52.