Most analytical and semi-analytical models for pumping-induced land subsidence invoke the simplifying assumptions regarding characteristics of geomaterials, as well as the pattern of drawdown response to pumping. This paper presents an analytical solution for one-dimensional consolidation of the multilayered soil due to groundwater drawdown, in which viscoelastic property and time-dependent drawdown are taken into account. The presented solution is developed by using the boundary transformation techniques. The validity of the proposed solution is verified by comparing with a degenerated case for a single layer, as well as with the numerical solutions and experimental results for a two-layer system. The difference between the average consolidation degree Up defined by hydraulic head and that Us defined by total settlement is discussed. The detailed parametric studies are conducted to reveal the effects of viscoelastic properties and drawdown patterns on the consolidation process. It is revealed that while the effect of different drawdown response patterns is significant during the early-intermediate stages of consolidation, the viscoelastic properties may have a more dominant influence on long-term consolidation behavior, depending on the values of the material parameters, which are reflected in both the deformation process of soil layers and the dissipation of excess pore-water pressure.