In this paper, a novel two-dimensional (2D) non-stationary wideband geometry-based stochastic model (GBSM) for massive multiple-input multiple-output (MIMO) communication systems is proposed. Key characteristics of massive MIMO channels such as near field effects and cluster evolution along the array are addressed in this model. Near field effects are modelled by a second-order approximation to spherical wavefronts, i.e., parabolic wavefronts, leading to linear drifts of the angles of multipath components (MPCs) and non-stationarity along the array. Cluster evolution along the array involving cluster (dis)appearance and smooth average power variations is considered. Cluster (dis)appearance is modeled by a two-state Markov process and smooth average power variations are modelled by a spatial lognormal process. Statistical properties of the channel model such as time autocorrelation function (ACF), spatial cross-correlation function (CCF), and cluster average power and Rician factor variations over the array are derived. Finally, simulation results are presented and analyzed, demonstrating that parabolic wavefronts and cluster soft evolution are good candidates to model important massive MIMO channel characteristics.