Three deep sea Δ18O records are compared in the temporal and spectral domain with each other, a simple ice sheet and bedrock climate model simulation and insolation anomalies at 65 ¿N. The model is one-dimensional; it represents the accumulation and ablation processes of ice, its viscoplastic flow, and the isostatic adjustment of bedrock under it, being forced solely by orbital insolation anomalies. The objective of the comparison is to determine the importance of internal variability arising from feedback between the ice sheets and bedrock through the late Pliocene and Pleistocene in the simulations and in the climate record. It is found that a considerable part of the spectral variability at suborbital frequencies, i.e., below the precessional and obliquity frequencies, as seen in the sediment cores, is simulated by the model. Simulated variability near 100 kyr has its period determined by the difference between two precessional frequencies. The model's internal variability mechanism has a period near 30--50 kyr but is highly damped. Still, it enhances suborbital variability near 100 kyr and 300 kyr and appears to induce nonstationarity of the orbital and suborbital peaks on even longer time scales. Since the insolation forcing exhibits no noticeable peaks at suborbital or lower frequencies, it appears that the nonstationarity in the model simulations as well as the suborbital peaks themselves are due to coupling between model nonlinearities and the orbital forcing. Since the variability at 100 kyr in the deep sea records is correlated with the simulated variability, such a control by the orbital forcing may be also operating in the climate system. The absence from the simulation of a trend toward larger ice sheets and colder climate seen throughout the Pleistocene in the sediment record is possibly due to slow tectonics and CO2 forcing changes not included in the model. ¿ American Geophysical Union 1993