I've updated the material properties to vary with temperature, mainly in the range of 90-140 K. Using the parametric sweep function to vary the input power of the heater, we get the eigenfreqencies' dependence on temperature to show up. The fractional dependence of 1.3e-5 /K around 123 K matches with what Aaron calculated in this elog entry, which is always a good sign that nothing is going horribly wrong. I've also added the flats to the silicon disc, for better accuracy. See the screenshots below showing the frequency shift with temperature.