We present the detailed study of the thermodynamics of vibrational modes in disordered elastic systems such as the Bragg glass phase of lattices pinned by quenched impurities. Our study and our results are valid within the (mean field) replica Gaussian variational method. We obtain an expression for the internal energy in the quantum regime as a function of the saddle point solution, which is then expanded in powers of ħ at low temperature T. In the calculation of the specific heat Cv a non trivial cancellation of the term linear in T occurs, explicitly checked to second order in ħ. The final result is Cv ∝T3 at low temperatures in dimension three and two. The prefactor is controlled by the pinning length. This result is discussed in connection with other analytical or numerical studies