For the period July 2003 to August 2010, the interplanetary coronal mass ejection (ICME) catalogue maintained by Richardson and Cane lists 106 Earth-directed events, which have been measured in-situ by plasma and field instruments onboard the ACE satellite. We present a statistical investigation of the Earths thermospheric neutral density response by means of accelerometer measurements collected by the GRACE satellites, which are available for 104 ICMEs in the data set, and its relation to various geomagnetic indices and characteristic ICME parameters such as the impact speed, southward magnetic field strength (Bz). The majority of ICMEs causes a distinct density enhancement in the thermosphere, with up to a factor of eight compared to the pre-event level. We find high correlations between ICME Bz and thermospheric density enhancements (~0.9), while the correlation with the ICME impact speed is somewhat smaller (~0.7). The geomagnetic indices revealing the highest correlations are Dst and SYM-H (~0.9), the lowest correlations are obtained for kp and AE (~0.7), which show a nonlinear relation with the thermospheric density enhancements. Separating the response for the shock sheath region and the magnetic structure of the ICME, we find that the Dst and SYM-H reveal a tighter relation to the Bz minimum in the magnetic structure of the ICME, whereas the polar cap indices show higher correlations with the Bz minimum in the shock sheath region. Since the strength of the Bz component - either in the sheath or the magnetic structure of the ICME - is highly correlated (~0.9) with the neutral density enhancement, we discuss the possibility of satellite orbital decay estimates based on magnetic field measurements at L1, i.e. before the ICME hits the Earths magnetosphere. This will further stimulate progress in space weather understanding and applications regarding satellite operations.