Turbulence in the upper ocean generates fluctuations in temperature and salinity, which result in variations in inherent optical properties (IOPs) and further change the underwater light field. A simulation-based study is performed for the radiative transfer (RT) of natural light in the turbulent flows in the upper ocean. For a canonical problem of turbulent shear flow interacting with the sea surface with and without surface waves, large-eddy simulations are performed for fluid motions and the transport of temperature and salinity. Based on the resolved turbulence temperature and salinity fields, IOP variations are quantified, and the inhomogeneous RT equation is then simulated using a Monte Carlo method. Through the simulations of a variety of cases with different flow, temperature, and salinity conditions, the statistics of downwelling irradiance are quantified and analyzed. It is found that the vertical profile of the mean downwelling irradiance is mainly determined by the vertical structure of the mean values of the IOPs; and turbulence effect is manifested in the horizontal variations of the downwelling irradiance. The magnitude of the irradiance variation is governed by the differences in the temperature and salinity between their values at the surface and in the deep region. In the presence of surface waves, the irradiance variation is enhanced due to the surface deformation, which is also largely affected by wave-turbulence interaction. The LES and inhomogeneous RT simulation may provide a useful tool for the characterization of upper-ocean turbulence processes based on underwater RT measurements.