Mean values and differential distributions of event-shape variables have been studied in neutral current deep inelastic scattering using an integrated {luminosity} of 82.2 pb$^{-1}$ collected with the ZEUS detector at HERA. The kinematic range was $80 < Q^2 < 20 480\gev^2$ and $0.0024 < x < 0.6$, where $Q^2$ is the virtuality of the exchanged boson and $x$ is the Bjorken variable. The data are compared with a model based on a combination of next-to-leading-order QCD calculations with next-to-leading-logarithm corrections and the Dokshitzer-Webber non-perturbative power corrections. The power-correction method provides a reasonable description of the data for all event-shape variables studied. Nevertheless, the lack of consistency of the determination of $��_s$ and of the non-perturbative parameter of the model, $\albar$, suggests the importance of higher-order processes that are not yet included in the model.

35 pages, 11 figures