Volcano deformation has often been analysed in terms of crustal elasticity, where the temporal volume change of the deformation source is required to explain any transient behaviour of ground displacement. However, the response of viscoelastic crust, both to magma discharge and eventual subsequent recharge, is not negligible because the effective crustal viscosity is expected to be lowered by higher geothermal gradient at volcanoes. The present study, using a 3-D finite element model with a linear elastic layer underlain by a linear Maxwell viscoelastic layer with spatially uniform viscosity, examines the viscoelastic crustal response to both magma discharge and recharge into the upper crust. The model considers a sill-like magma chamber emplaced at the bottom of the elastic layer inflating at a constant rate, during which the magma discharge due to an eruption occurs instantaneously all at once. Numerical experiments have revealed that the viscoelasticity in response to magma discharge promotes surface uplift. The uplift depends on how long the pre-eruptive magma charge lasted, that is, how much viscoelastic relaxation in response to the magma charge occurred in pre-eruption period. If the relaxation was significant, the uplift due to viscoelastic relaxation is greater than that due to magma recharge earlier in post-eruption period. On the contrary, at later times, the uplift due to magma recharge exceeds the relaxation-derived one. However, since the viscoelastic response to magma recharge detracts the uplift, the post-eruptive magma volume is recovered before the surface elevation is recovered. This viscoelastic model behaviour is then adopted for the geodetic data in the Aira caldera, southern Kyushu, Japan, giving important constraints on the model parameters. A model with constant magma supply rate explains the observed exponential-like recovery of the surface elevation after the 1914 eruption if the following conditions are satisfied: (1) the effective crustal viscosity is ~5x10^18 Pa s, and (2) magma emplacement at a depth of ~10-15 km is initiated ~50 years before the eruption, with a constant supply rate of at least ~0.01-0.02 km^3/yr. The constrained model predicts that the accumulated volume of magma during the first 50 years in the post-eruption period is lower than that estimated by elastic model, but is on the contrary larger thereafter. In addition, the surface uplift immediately before the 1914 eruption is predicted to be ~3 times lower than that supposed in previous studies. Our results suggest that geodetic signals observed before and after the magma discharge due to an eruption need to be revisited with respect to the crustal viscoelasticity.