Landsat has a history of use in the diagnosis of land surface phenology, vegetation disturbance, and their impacts on numerous forest biological processes. Studies have connected remote sensing-based phenology to surface climatological patterns, often using average temperatures and derived growing degree day accumulations. I present a detailed examination of remotely sensed forest phenology in the region of western Lake Superior, USA, based on a comprehensive climatological assessment and 1984-2013 Landsat imagery. I use this climatology to explain both the mean annual land surface phenological cycle and its interannual variability in temperate mixed forests. I assess long-term climatological means, trends, and interannual variability for the study period using available weather station data, focusing on numerous basic and derived climate indicators: seasonal and annual temperature and precipitation, the traditionally defined frost-free growing season, and a newly defined metric of the climatological growing season: the warm-season plateau in accumulated chilling days. Results indicate +0.56°C regional warming during the 30-year study period, with cooler springs (–1.26°C) and significant autumn warming (+1.54°C). The duration of the climatological growing season has increased +0.27 days/y, extending primarily into autumn. Summer precipitation in my study area declined by an average –0.34 cm/y, potentially leading to moisture stress that can impair vegetation carbon uptake rates and can render the forest more vulnerable to disturbance. Many changes in temperature, precipitation, and climatological growing season are most prominent in locations where Lake Superior exerts a strong hydroclimatological influence, especially the Minnesota shoreline and in forest areas downwind (southeast) of the lake. I then develop and demonstrate a novel phenoclimatological modeling method, applied over five Landsat footprints across my study area, that combines (1) diagnosis of the mean phenological cycle from remote sensing observations with (2) a partial-least-squares regression (PLSR) approach to modeling vegetation index residuals using these numerous meteorological and climatological observations. While the mean phenology often used to inform land surface models in meteorological and climate modeling systems may explain 50-70% of the observed phenological variability, this mixed modeling approach can explain more than 90% of the variability in phenological observations based on long-term Landsat records for this region.