The accurate retrieval of land surface vegetation properties under varying environmental conditions from time series of moderately high spatial resolution satellite observations is challenging. By coupling various Radiative Transfer (RT) models one can describe the soil, vegetation and atmosphere contributions in a “bottom-up” approach and, thereby, simulate top-of-atmosphere (TOA) spectral radiance data comparable to satellite-observed TOA radiances. This makes it possible to retrieve vegetation properties directly from TOA radiances rather than from atmospherically corrected top-of-canopy (TOC) reflectance data. The advantages of this approach are that a separate atmospheric correction of the satellite images is not necessary, and that the anisotropic surface reflection can also be taken into account effectively. In this study, we coupled various RT models, including the brightness – shape – moisture (BSM) reflectance model of the soil, the optical radiative transfer (RTMo) model of vegetation and the ‘MODerate resolution atmospheric TRANsmission’ (MODTRAN) model of the atmosphere, to simulate an annual time series of Landsat satellite TOA radiances observed during a drought episode in California Mediterranean grasslands in 2004. The inversion of this coupled system through an optimization technique against Landsat TOA radiances resulted in direct retrieval of vegetation properties. We accommodated the surface anisotropic reflection in our coupled modeling and also defined a novel anisotropy index to quantitatively express the importance of this phenomenon in satellite image analysis for the first time. Our study showed that the coupled use of RT models was able to accurately reproduce the time series of observed TOA radiances collected under varying soil moisture contents during an extended drought episode. The proposed coupling approach is useful for successful retrieval of vegetation properties from time series of satellite TOA radiance data to produce maps of land surface properties and to monitor vegetation properties variations in an operational straightforward way. The approach can also be easily adapted for conducting multi-sensor time series studies, creating a much denser temporal sampling than would be possible for separate single sensors.