  We perform three dimensional radiation hydrodynamic simulations of the
structure and dynamics of radiation dominated envelopes of massive stars at the
location of the iron opacity peak. One dimensional hydrostatic calculations
predict an unstable density inversion at this location, whereas our simulations
reveal a complex interplay of convective and radiative transport whose behavior
depends on the ratio of the photon diffusion time to the dynamical time. The
latter is set by the ratio of the optical depth per pressure scale height,
$\tau_0$, to $\tau_c=c/c_g$, where $c_g \approx$ 50 km/s is the isothermal
sound speed in the gas alone. When $\tau_0 \gg \tau_c$, convection reduces the
radiation acceleration and removes the density inversion. The turbulent energy
transport in the simulations agrees with mixing length theory and provides its
first numerical calibration in the radiation dominated regime. When $\tau_0 \ll
\tau_c$, convection becomes inefficient and the turbulent energy transport is
negligible. The turbulent velocities exceed $c_g$, driving shocks and large
density fluctuations that allow photons to preferentially diffuse out through
low-density regions. However, the effective radiation acceleration is still
larger than the gravitational acceleration so that the time average density
profile contains a modest density inversion. In addition, the simulated
envelope undergoes large-scale oscillations with periods of a few hours. The
turbulent velocity field may affect the broadening of spectral lines and
therefore stellar rotation measurements in massive stars, while the time
variable outer atmosphere could lead to variations in their mass loss and
stellar radius.
