We present a new set of horizontal branch (HB) models computed with the MESA stellar evolution code. The models adopt alpha-enhanced Asplund et al. metal mixtures and include the gravitational settling of He. They are used in our HB population synthesis tool to generate theoretical distributions of HB stars in order to describe the multiple stellar populations in the globular clusters 47 Tuc, M3, and M13. The observed HB in 47 Tuc is reproduced very well by our simulations for [Fe/H] = -0.70 and [alpha/Fe] = +0.4 if the initial helium mass fraction varies by Delta Y-0 similar to 0.03, and approximately 21%, 37%, and 42% of the stars have Y-0 = 0.257, 0.270, and 0.287, respectively. These simulations yield (m - M)(V) = 13.27, implying an age near 13.0 Gyr. In the case of M3 and M13, our synthetic HBs for [Fe/H] = -1.55 and [alpha/Fe] = 0.4 match the observed ones quite well if M3 has Delta Y-0 similar to 0.01 and (m - M)(V) = 15.02, resulting in an age of 12.6 Gyr, whereas M13 has Delta Y-0 similar to 0.08 and (m - M)(V) = 14.42, implying an age of 12.9 Gyr. Mass loss during giant branch evolution and Delta Y-0 appear to be the primary second parameters for M3 and M13. New observations for seven of the nine known RR Lyrae in M13 are also reported. Surprisingly, periods predicted for the c-type variables tend to be too high (by up to similar to 0.1 days).