We describe and validate a Monte Carlo model to track photons over the full range of solar wavelengths as they travel into optically thick Antarctic blue ice. The model considers both reflection and transmission of radiation at the surface of blue ice, scattering by air bubbles within it and spectral absorption due to the ice. The ice surface is treated as planar whilst bubbles are considered as spherical scattering centres using the Henyey-Greenstein approximation. Using bubble radii and number concentrations that are representative of Antarctic blue ice, we calculate spectral albedos and spectrally-integrated downwelling and upwelling radiative fluxes as functions of depth and find that, relative to the incident irradiance, there is a marked subsurface enhancement in the downwelling flux and accordingly also in the mean irradiance. This is due to the interaction between the refractive air-ice interface and the scattering interior and is particularly notable at blue and UV wavelengths which correspond to the minimum of the absorption spectrum of ice. In contrast the absorption path length at IR wavelengths is short and consequently the attenuation is more complex than can be described by a simple Lambert-Beer style exponential decay law — instead we present a triple exponential fit to the net irradiance against depth. We find that there is a moderate dependence on the solar zenith angle and surface conditions such as altitude and cloud optical depth. Representative broadband albedos for blue ice are calculated in the range 0.585 to 0.621. For macroscopic absorbing inclusions we observe both geometry- and size-dependent self-shadowing that reduces the fractional irradiance incident on an inclusion’s surface. Despite this, the inclusions act as local photon sinks and are subject to fluxes that are several times the magnitude of the single scattering contribution. Such enhancement may have consequences for the energy budget in regions of the cryosphere where particulates are present near the surface. These results also have particular relevance to measurements of the internal radiation field: account must be taken of both self-shadowing and the optical effect of introducing the detector. Turning to the particular example of englacial meteorites, our modelling predicts iron meteorites to reside at much reduced depths than previously suggested in the literature (< 10 cm vs ~40 cm), and further shows a size dependency that may explain the observed bias in their Antarctic size distribution.