This work is concerned with asymptotic approximations of the semi-classical Schr\"odinger equation in periodic media using Gaussian beams. For the underlying equation, subject to a highly oscillatory initial data, a hybrid of the Gaussian beam approximation and homogenization leads to the Bloch eigenvalue problem and associated evolution equations for Gaussian beam components in each Bloch band. We formulate a superposition of Bloch-band based Gaussian beams to generate high frequency approximate solutions to the original wave field. For initial data of a sum of finite number of band eigen-functions, we prove that the first-order Gaussian beam superposition converges to the original wave field at a rate of ϵ1/2, with ϵ the semiclassically scaled constant, as long as the initial data for Gaussian beam components in each band are prepared with same order of error or smaller. For a natural choice of initial approximation, a rate of ϵ1/2 of initial error is verified.