Hierarchical equations of motion method applied to nonequilibrium heat transport in model molecular junctions: Transient heat current and high-order moments of the current operator