We present radiation transfer (RT) simulations of evolutionary sequences of massive protostars forming from massive dense cores in environments of high surface densities. The protostellar evolution is calculated with a detailed multi-zone model, with
the accretion rate regulated by feedback from an evolving disk-wind outflow cavity. Disk and envelope evolutions are calculated self-consistently. In this framework, an evolutionary track is determined by three environmental initial conditions: the initial core mass M_c, the mean surface density of the ambient star-forming clump Sigma_cl, and the rotational-to-gravitational energy ratio of the initial core, beta_c. Evolutionary sequences with various M_c, Sigma_cl, beta_c are constructed. We find that in a fiducial model with M_c=60Msun, Sigma_cl=1 g/cm^2 and beta_c=0.02, the final star formation efficiency >~0.43. For each evolutionary track, RT simulations are performed at selected stages, with temperature profiles, SEDs, and images produced. At a given stage the envelope temperature is highly dependent on Sigma_cl, but only weakly dependent on M_c. The SED and MIR images depend sensitively on the evolving outflow cavity, which gradually wides as the protostar grows. The fluxes at <~100 microns increase dramatically, and the far-IR peaks move to shorter wavelengths. We find that, despite scatter caused by different M_c, Sigma_cl, beta, and inclinations, sources at a given evolutionary stage appear in similar regions on color-color diagrams, especially when using colors at >~ 70 microns, where the scatter due to the inclination is minimized, implying that such diagrams can be useful diagnostic tools of evolutionary stages of massive protostars. We discuss how intensity profiles along or perpendicular to the outflow axis are affected by environmental conditions and source evolution.
