The S shaped log log drawdown time curve typical of pumping tests in unconfined aquifers is reinvestigated via numerical experiments. We examine the temporal and spatial evolution of the rate of change in storage in an unconfined aquifer during pumping. This evolution is related to the transition of water release mechanisms from the expansion of water and compaction of the porous medium to the drainage of water from the unsaturated zone above the initial water table and initially saturated pores as the water table falls during the pumping of the aquifer. We conclude that the transition of the water release mechanisms and vertical flow are the cause of the S shaped drawdown time. Cross-correlation analysis is then employed to examine the relationship between hydraulic properties of an unconfined aquifer and pressure observations. The analysis reveals that head observed in the saturated zone at late times along a streamline is positively correlated with the conductivity (K(s)) of the region upstream of the observation location, and negatively correlated with the K(s) of the region downstream of the observation location along the same streamline. Besides, head observations in the saturated zone at the early time are positively correlated with specific storage (S(s)) in a narrow region between the observation and pumping locations. At intermediate and late times, the head positively correlates with the heterogeneity of α (pore-size distribution parameter) in a thin disk-shaped unsaturated region above the pumping and observation locations. Saturated water content θ(s) in the vadose zone directly above the pumping and monitoring locations is found positively correlated with the head observations during the intermediate times and late times.In the end, a stochastic inverse estimation is conducted to jointly interpret a sequential pumping test in a three dimensional unconfined aquifer. K(s), S(s), θ(s) and α are estimated at the same time. The estimated results capture the pattern of the heterogeneous parameters as well as the details with a smooth distribution. The estimated heterogeneous parameter fields produce better head predictions than the traditional homogeneous method.