In this paper, we investigate a sparse channel estimation problem for broadband massive multiple-input-multiple-output orthogonal frequency division multiplexing (MIMO-OFDM) systems. We propose a hidden Markov model to capture the structured sparsity and temporal dependency characteristic of massive MIMO-OFDM channels in the angle-delay domain, and this probability model exhibits extensive adaptability to different realistic propagation scenarios. Then we solve the channel estimation problem based on a novel optimization framework named constrained Bethe free energy (BFE) minimization, which is valid for a generic statistical model. Under this systematic theoretical framework, a hierarchical hybrid message passing (HHMP) algorithm is proposed to track dynamic channel parameters recursively. The proposed method can adaptively learn the sparse structure and temporal correlation of multiuser channels without requiring the knowledge of hidden Markov channel parameters. Numerical simulations demonstrate that the proposed HHMP algorithm can accurately estimate angle-delay domain channels with reduced iteration times and pilot overhead.